diff --git a/README.md b/README.md index e69de29..6bd0fc7 100644 --- a/README.md +++ b/README.md @@ -0,0 +1,24 @@ +# Project 1: Regression and Resampling + +This project was authored by **Lars Bogner**. + +## Project Structure + +- **`src/`**: Contains the main code for the project. The majority of the implementation is contained in the `pyoptim` library. + - **`main.ipynb`**: Jupyter Notebook used for figure generation. +- **`report/`**: Contains the LaTeX source files for the project report. + +## Installation + +To install all dependencies, use the following command: + +```bash +uv install +``` +This uses the `uv` tool to manage the environment and dependencies. Refer to the [uv documentation](https://uv.readthedocs.io/en/latest/) for more details. + +## Usage + +1. Navigate to the `src/` directory. +2. Open and run `main.ipynb` to generate figures. +3. The generated figures will be saved in the `figures/` directory. \ No newline at end of file diff --git a/report/chapters/methods.tex b/report/chapters/methods.tex index f519b2d..0029875 100644 --- a/report/chapters/methods.tex +++ b/report/chapters/methods.tex @@ -57,7 +57,7 @@ where the $L_1$ norm is defined as $\norm{\vec \theta}_1 = \sum_{i=0}^{p-1} |\th \begin{equation} \label{eq:ridge_dof} n_\text{effective} = \sum_{i=1}^{n_\text{features}} \frac{d_i^2}{d_i^2 + \lambda}, \end{equation} -where $d_i$ are the singular values of the design matrix $X$. Thus regularization effectively reduces the number of parameters in the model, leading to a bias-variance tradeoff that can be tuned by the regularization parameter $\lambda$. While OLS models can be analytically optimized using \cref{eq:ols_solution} for the case where $X^TX$ is invertible, the same is true for Ridge regression in the case where $X^TX + \lambda I$ is invertible. The optimal parameters for Ridge regression can be found by solving the modified normal equations\cite{elstner_lecture_2025} +where $d_i$ are the singular values of the design matrix $X$. Thus, regularization effectively reduces the number of parameters in the model, leading to a bias-variance tradeoff that can be tuned by the regularization parameter $\lambda$. While OLS models can be analytically optimized using \cref{eq:ols_solution} for the case where $X^TX$ is invertible, the same is true for Ridge regression in the case where $X^TX + \lambda I$ is invertible. The optimal parameters for Ridge regression can be found by solving the modified normal equations\cite{elstner_lecture_2025} \begin{equation} \label{eq:ridge_solution} \vec \theta_\mathrm{Ridge} = (X^TX + \lambda I)^{-1}X^T \vec y. \end{equation} @@ -179,7 +179,7 @@ The process of k-fold cross-validation is illustrated in \cref{fig:crossvalidati (6.8*\squaresize,\ybottom) node[midway,xshift=1.5cm]{$\displaystyle \frac{1}{5}\sum_{i=1}^{5}\mathrm{MSE}_i$}; \end{tikzpicture} - \caption{Illustration of 5-fold cross-validation. Each row corresponds to one fold: the turquoise block marks the held-out test set while the white blocks form the training set. The models error $\mathrm{MSE}_i$ is computed for each fold, and the final performance is the average of all $\mathrm{MSE}_i$.} + \caption{Illustration of 5-fold cross-validation. Each row corresponds to one fold: the turquoise block marks the held-out test set while the white blocks form the training set. The models' error $\mathrm{MSE}_i$ is computed for each fold, and the final performance is the average of all $\mathrm{MSE}_i$.} \label{fig:crossvalidation} \end{figure} @@ -189,7 +189,7 @@ The process of k-fold cross-validation is illustrated in \cref{fig:crossvalidati The methods described in the previous section have been implemented in Python. The implementation is structured in a modular way, allowing to easily switch between different optimization algorithms and resampling methods. The code is available on GitHub at \url{https://github.uio.no/larsbog/FYSSTK-Project1}. The following libraries have been used in the implementation: \begin{description} - \item[\texttt{numpy}\cite{harris_array_2020}] It is a fundamental package for scientific computing in Python. It provides support for large, multi-dimensional arrays and matrices, along with a collection of mathematical functions to operate on these arrays efficiently. It is used in this project for all matrix and vector operations, as well as for generating random numbers and performing statistical calculations. + \item[\texttt{numpy}\cite{harris_array_2020}] It is a fundamental package for scientific computing in Python. It provides support for large, multidimensional arrays and matrices, along with a collection of mathematical functions to operate on these arrays efficiently. It is used in this project for all matrix and vector operations, as well as for generating random numbers and performing statistical calculations. \item[\texttt{scikit-learn}\cite{pedregosa_scikit-learn_2011}] It is a machine learning library for Python that provides simple and efficient tools for data mining and data analysis. It includes implementations of various machine learning techniques, like data manipulation and performance evaluation metrics. It is used in this project for rescaling the data, splitting the data into training and test sets, and for implementing the resampling methods. Furthermore, it provides the implementation of the mean squared error (MSE) metric used to evaluate the performance of the models. \item[\texttt{matplotlib}\cite{hunter_matplotlib_2007}] It is a plotting library for Python that provides a wide range of tools for creating different types of plots and visualizations. It is used in this project to visualize the results of the different methods and to create plots for the report. @@ -206,6 +206,6 @@ The Runge function is defined as \begin{equation} \label{eq:runge_function} f(x) = \frac{1}{1 + 25x^2}, \quad x \in [-1, 1]. \end{equation} -The $x$-values are sampled uniformly in the interval $[-1, 1]$, and the corresponding $y$-values are generated by adding Gaussian noise with mean 0 and standard deviation $\sigma = 1$, unless otherwise specified. To evaluate the performance in dependance of the degrees of freedom, the $x$-values are expanded into multiple polynomial features up to a specified degree: $X_i = \left\{x_i^0, x_i^1, \dots, x_i^p\right\}$. All features and outputs are then rescaled to have zero mean and unit variance using the \texttt{StandardScaler} from \texttt{scikit-learn}. The mean and variance of the training data are used for rescaling the test data to avoid data leakage. -There are analytical implementations of the OLS and Ridge regression methods, as described in \cref{eq:ols_solution} and \cref{eq:ridge_solution}. Furthermore numerical optimizers are implemented from scratch based on a class inheritance structure to reduce code duplication. The base class \texttt{GradientDescent} implements a fitting procedure consisting of a precomputation step, an optimization procedure over multiple iterations and a memory of the cost function values. The different optimization algorithms, i.e. vanilla gradient descent, momentum based gradient descent, Adagrad, RMSProp and Adam, are implemented as subclasses that override the update rule for the parameters $\vec \theta$ in each iteration. The cost functions for OLS, Ridge and Lasso regression are also implemented as separate classes that provide methods to calculate the cost and its gradient. The stochastic gradient descent methods are implemented as further subclasses that modify the data used for fitting in each iteration. In the stochastic gradient descent method, the cost history property is furthermore modified to return the per epoch average of the cost function values. The resampling methods, i.e. bootstrapping and k-fold cross-validation, are implemented as methods with identical interfaces which provide multiple sets of training and test data. This allows to easily switch between the different resampling methods when evaluating the performance of a model. +The $x$-values are sampled uniformly in the interval $[-1, 1]$, and the corresponding $y$-values are generated by adding Gaussian noise with mean 0 and standard deviation $\sigma = 1$, unless otherwise specified. To evaluate the performance in dependence of the degrees of freedom, the $x$-values are expanded into multiple polynomial features up to a specified degree: $X_i = \left\{x_i^0, x_i^1, \dots, x_i^p\right\}$. All features and outputs are then rescaled to have zero mean and unit variance using the \texttt{StandardScaler} from \texttt{scikit-learn}. The mean and variance of the training data are used for rescaling the test data to avoid data leakage. +There are analytical implementations of the OLS and Ridge regression methods, as described in \cref{eq:ols_solution} and \cref{eq:ridge_solution}. Furthermore, numerical optimizers are implemented from scratch based on a class inheritance structure to reduce code duplication. The base class \texttt{GradientDescent} implements a fitting procedure consisting of a precomputation step, an optimization procedure over multiple iterations and a memory of the cost function values. The different optimization algorithms, i.e. vanilla gradient descent, momentum based gradient descent, Adagrad, RMSProp and Adam, are implemented as subclasses that override the update rule for the parameters $\vec \theta$ in each iteration. The cost functions for OLS, Ridge and Lasso regression are also implemented as separate classes that provide methods to calculate the cost and its gradient. The stochastic gradient descent methods are implemented as further subclasses that modify the data used for fitting in each iteration. In the stochastic gradient descent method, the cost history property is furthermore modified to return the per epoch average of the cost function values. The resampling methods, i.e. bootstrapping and k-fold cross-validation, are implemented as methods with identical interfaces which provide multiple sets of training and test data. This allows to easily switch between the different resampling methods when evaluating the performance of a model. diff --git a/report/chapters/results.tex b/report/chapters/results.tex index 1c3502b..b60eaaa 100644 --- a/report/chapters/results.tex +++ b/report/chapters/results.tex @@ -20,7 +20,7 @@ How the mean squared error (MSE) and the $R^2$ score depend on the polynomial de \label{fig:ols_params} \end{figure} -In the limit of small datasets this issue becomes even more relevant. As the number of parameters approaches the number of data points, the model can fit the data perfectly, leading to a MSE of zero. However, this is usually not a desirable outcome, as the model will not generalize well to new data, as it has simply \textquote{memorized} the noise in the training data. This is a classic example of overfitting, where the model performs well on the training data but poorly on unseen data. In such cases, it is crucial to use techniques to prevent overfitting, such as regularization or cross-validation. +In the limit of small datasets this issue becomes even more relevant. As the number of parameters approaches the number of data points, the model can fit the data perfectly, leading to an MSE of zero. However, this is usually not a desirable outcome, as the model will not generalize well to new data, as it has simply \textquote{memorized} the noise in the training data. This is a classic example of overfitting, where the model performs well on the training data but poorly on unseen data. In such cases, it is crucial to use techniques to prevent overfitting, such as regularization or cross-validation. \subsubsection{Regularization Techniques} \label{sec:results_reg_techniques} @@ -59,7 +59,7 @@ With an increasing regularization strength $\lambda$ the parameter values are pu -However a second aspect of SGD is visible in \cref{fig:optimization_performance}. While GD converges smoothly to the optimal solution, SGD shows significant oscillations around the optimal solution. This is a direct consequence of the increased variance in the gradient estimate, as the gradient is only computed over a small subset of the data. This can lead to slower convergence and even no convergence at all, if the batch size is too small or the noise in the data is too high. The intrinsic noise in the gradient estimate leads also to the issue, that SGD may be unable to reach an optimal solution with a large amount of noise in the training data. To counteract this problem, the noise in the dataset has been reduced by one order of magnitude compared to the previous sections. However, even with a reduced noise level, SGD is unable to reach the optimal solution. Further fine tuning of the new hyperparameters introduced by SGD, i.e. the batch size and the batches per epoch, could lead to better results. However, this is a non-trivial task and is highly dependent on the problem at hand. To arrive at a robust solution it might be necessary to increase the batch size until reaching the limit of full-batch GD, which would defeat the purpose of using SGD in the first place. +However, a second aspect of SGD is visible in \cref{fig:optimization_performance}. While GD converges smoothly to the optimal solution, SGD shows significant oscillations around the optimal solution. This is a direct consequence of the increased variance in the gradient estimate, as the gradient is only computed over a small subset of the data. This can lead to slower convergence and even no convergence at all, if the batch size is too small or the noise in the data is too high. The intrinsic noise in the gradient estimate leads also to the issue, that SGD may be unable to reach an optimal solution with a large amount of noise in the training data. To counteract this problem, the noise in the dataset has been reduced by one order of magnitude compared to the previous sections. However, even with a reduced noise level, SGD is unable to reach the optimal solution. Further, fine-tuning of the new hyperparameters introduced by SGD, i.e. the batch size and the batches per epoch, could lead to better results. However, this is a non-trivial task and is highly dependent on the problem at hand. To arrive at a robust solution it might be necessary to increase the batch size until reaching the limit of full-batch GD, which would defeat the purpose of using SGD in the first place. \subsection{Effect of Learning Rate on Convergence} \label{sec:results_learning_rate} @@ -77,7 +77,7 @@ In \cref{fig:learning_rates} it is easily observable, that the convergence towar \subsection{Alternative Optimization Algorithms} \label{sec:results_alt_opt_algos} -A different approach to limit the computational cost of fitting models with a high number of parameters is to use alternative optimization algorithms that reduce the number of gradient evaluations before reaching convergence. The optimized algorithms as introduced in \Cref{sec:numerical_optimization} are evaluated in regards to their convergence behavior on the same dataset and model. +A different approach to limit the computational cost of fitting models with a high number of parameters is to use alternative optimization algorithms that reduce the number of gradient evaluations before reaching convergence. The optimized algorithms as introduced in \Cref{sec:numerical_optimization} are evaluated with regard to their convergence behavior on the same dataset and model. \begin{figure} \centering @@ -93,13 +93,13 @@ A different approach to limit the computational cost of fitting models with a hi \begin{figure} \centering \includegraphics[width=\columnwidth]{../figures/bias_variance_tradeoff_combined_plot.pdf} - \caption{Bias-variance decomposition of the mean squared error (MSE) for polynomial fits of different degrees to noisy data from the Runge function using OLS. The Bias-Variance decomposition is calculated from bootstraped samples. The MSE is also calculated using 5-fold cross-validation and bootstrap resampling techniques. $\N=\num{80}, \Nt = \num{20}$.} + \caption{Bias-variance decomposition of the mean squared error (MSE) for polynomial fits of different degrees to noisy data from the Runge function using OLS. The Bias-Variance decomposition is calculated from bootstrapped samples. The MSE is also calculated using 5-fold cross-validation and bootstrap resampling techniques. $\N=\num{80}, \Nt = \num{20}$.} \label{fig:bias_variance_tradeoff} \end{figure} To evaluate the bias-variance tradeoff of the optimized models, the MSE is decomposed into its bias and variance components as introduced in \cref{sec:bias_variance_tradeoff}. The results are shown in \cref{fig:bias_variance_tradeoff}. To accurately depict the bias-variance tradeoff, there were two sets of test data used. To compute the bias, no noise was added to the test data, while for the rest of the computations the same noise level as in the training data was used. It is visible that for increasing polynomial degrees the variance increases, while the bias stays almost constant. This is a direct consequence of the increased flexibility of the model with increasing polynomial degree. To ensure numerical stability, Ridge regression with a very small regularization strength of $\lambda = \num{1e-10}$ has been used to compute the bias-variance decomposition in an approximation of OLS. This ensures that the matrix inversion in the computation of the model parameters is numerically stable, while the regularization term has almost no effect on the model performance. The optimal model complexity is usually found at the point where the sum of bias and variance is minimal. Usually this coincides with the point of the minimum of the mean squared error. How the MSE behaves on average for $k$-fold cross-validation and bootstrap resampling techniques is also shown in \cref{fig:bias_variance_tradeoff}. -Both resampling techniques show a minimum in the MSE for for the polynomial degree of 2. The crossvalidation technique shows a much more constant MSE value over the entire range of polynomial degrees, while the bootstrap technique shows a more pronounced minimum. This indicates that the bootstrap technique is more sensitive to the choice of polynomial degree and thus provides a better estimate in the case of hyperparameter tuning at the cost of a much higher computational cost. +Both resampling techniques show a minimum in the MSE for the polynomial degree of 2. The crossvalidation technique shows a much more constant MSE value over the entire range of polynomial degrees, while the bootstrap technique shows a more pronounced minimum. This indicates that the bootstrap technique is more sensitive to the choice of polynomial degree and thus provides a better estimate in the case of hyperparameter tuning at the cost of a much higher computational cost. This trend of optimal model complexity being achieved somewhere between the high bias case of low parameter count and the high variance case of high parameter count is a general property of machine learning models and is not limited to polynomial regression with OLS cost functions. As seen in \cref{fig:kfold_mse_comp} this trend of a minimum in the MSE is also observable for Ridge and Lasso regression. diff --git a/report/ex39.pdf b/report/ex39.pdf new file mode 100644 index 0000000..c474df7 Binary files /dev/null and b/report/ex39.pdf differ diff --git a/report/main.pdf b/report/main.pdf index 3759407..69275a3 100644 Binary files a/report/main.pdf and b/report/main.pdf differ diff --git a/report/main.tex b/report/main.tex index 9222473..01e23f7 100644 --- a/report/main.tex +++ b/report/main.tex @@ -1,5 +1,5 @@ -\documentclass[aps,rmp,reprint,amsmath,amssymb,graphicx,longbibliography]{revtex4-1} +\documentclass[aps,rmp,reprint,amsmath,amssymb,graphicx,longbibliography,twoside]{revtex4-1} \usepackage{bm} \usepackage{graphicx}