\documentclass[english,notitlepage,reprint,nofootinbib]{revtex4-2}

% ===== Imports =====
\usepackage[utf8]{inputenc}
\usepackage{physics,amssymb}
\usepackage{graphicx}
\usepackage{xcolor}
\usepackage{hyperref}
\usepackage{listings}
\usepackage{subfigure}
\usepackage{algorithm}
\usepackage[noend]{algpseudocode}
\usepackage{subfigure}
\usepackage{tikz}
\usepackage{float}

\hypersetup{
    colorlinks,
    linkcolor={red!50!black},
    citecolor={blue!50!black},
    urlcolor={blue!80!black}}

% ===========================================

\begin{document}

\title{Systematic Analysis of Feedforward Neural Network Architectures:\\Implementation, Optimization, and Comparative Performance Evaluation}
\author{Sebastian Koranda \quad Yoan Tutunarov \quad Vladislav Foss}
\date{\today}
\noaffiliation

\begin{abstract}
This study implements and analyzes feedforward neural networks (FFNNs) from first principles for both regression and classification tasks. We develop a flexible neural network framework using the backpropagation algorithm, supporting multiple activation functions (Sigmoid, ReLU, Leaky ReLU) and optimization methods (gradient descent, RMSprop, Adam). We compare against the performance of existing neural network libraries and validate the gradients against automatic differentiation.
\end{abstract}

\maketitle

\section{Introduction}
The primary motivation for this study lies in understanding the fundamental mechanics of neural networks through ground-up implementation rather than relying solely on high-level libraries. By developing our own FFNN code, we gain deep insights into the mathematical foundations of gradient-based optimization, activation function behavior, and the bias-variance tradeoff inherent in neural network architectures. This hands-on approach enables a more thorough understanding of how various hyperparameters—including learning rates, network depth, node counts, and regularization parameters—influence model performance and generalization capability.

Our investigation focuses on two complementary problem domains: regression using the challenging Runge function, and multiclass classification using the well-established MNIST digit recognition dataset. The Runge function serves as an excellent benchmark for regression analysis due to its highly nonlinear nature and tendency to create oscillatory behavior that challenges polynomial approximation methods. By comparing neural network performance against our previous OLS regression results from Project 1, we can quantitatively assess the advantages of nonlinear modeling approaches.

All code developed for this report is available via a GitHub repository.\footnote{\url{https://github.uio.no/sebasmk/fys-stk4155/tree/main/projects/project2}}

\subsection*{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{Methods}

This section describes the theoretical foundations and implementation details of our feedforward neural network framework. We present the mathematical formulation of the network architecture, the backpropagation algorithm for gradient computation, and the optimization strategies employed for training.

\subsection{Neural Network Architecture}

A feedforward neural network consists of an input layer, one or more hidden layers, and an output layer. Information flows unidirectionally from input to output without cycles or loops. For a network with $L$ layers (excluding the input layer), the forward propagation through layer $l$ is computed as:

\begin{equation}
    \mathbf{z}^{(l)} = \mathbf{W}^{(l)} \mathbf{a}^{(l-1)} + \mathbf{b}^{(l)}
\end{equation}
\begin{equation}
    \mathbf{a}^{(l)} = f^{(l)}(\mathbf{z}^{(l)})
\end{equation}

where $\mathbf{W}^{(l)}$ is the weight matrix connecting layer $l-1$ to layer $l$, $\mathbf{b}^{(l)}$ is the bias vector for layer $l$, $\mathbf{z}^{(l)}$ represents the pre-activation values, $\mathbf{a}^{(l)}$ denotes the activations (with $\mathbf{a}^{(0)} = \mathbf{X}$ being the input), and $f^{(l)}$ is the activation function for layer $l$.

The final output $\hat{\mathbf{y}} = \mathbf{a}^{(L)}$ represents the network's prediction. For regression tasks, we use a linear activation in the output layer, while for classification tasks, we employ softmax activation to produce probability distributions over classes.

\subsection{Weight Initialization}

Proper weight initialization is crucial for successful training of neural networks. We initialize weights using a Gaussian distribution with small standard deviation:

\begin{equation}
    W_{ij}^{(l)} \sim \mathcal{N}(0, 0.01^2)
\end{equation}

This small-scale initialization prevents saturation of activation functions and ensures gradients flow properly during backpropagation. Biases are initialized to zero:

\begin{equation}
    b_i^{(l)} = 0
\end{equation}

This initialization strategy is suitable for the sigmoid activation functions used in our hidden layers and helps avoid symmetry-breaking issues during training.

\subsection{Activation Functions}

Activation functions introduce non-linearity into the network, enabling it to learn complex patterns. We implemented and tested three different activation functions for hidden layers.

\subsubsection{Sigmoid Function}

The sigmoid function maps inputs to the range $(0, 1)$:

\begin{equation}
    \sigma(z) = \frac{1}{1 + e^{-z}}
\end{equation}

Its derivative, which is required for backpropagation, has the convenient form:

\begin{equation}
    \frac{d\sigma}{dz} = \sigma(z)(1 - \sigma(z))
\end{equation}

This allows efficient computation since the derivative can be expressed in terms of the function output itself. However, sigmoid activations can suffer from vanishing gradients for large $|z|$ values, as the derivative approaches zero.

\subsubsection{ReLU (Rectified Linear Unit)}

The ReLU activation is defined as:

\begin{equation}
    \text{ReLU}(z) = \max(0, z) = \begin{cases} 
    z & \text{if } z > 0 \\
    0 & \text{if } z \leq 0
    \end{cases}
\end{equation}

with derivative:

\begin{equation}
    \frac{d\text{ReLU}}{dz} = \begin{cases} 
    1 & \text{if } z > 0 \\
    0 & \text{if } z \leq 0
    \end{cases}
\end{equation}

ReLU addresses the vanishing gradient problem and is computationally efficient. However, neurons can "die" during training when they consistently output zero, preventing further learning.

\subsubsection{Leaky ReLU}

To address the dying ReLU problem, Leaky ReLU allows small negative values:

\begin{equation}
    \text{LeakyReLU}(z) = \begin{cases} 
    z & \text{if } z > 0 \\
    \alpha z & \text{if } z \leq 0
    \end{cases}
\end{equation}

with derivative:

\begin{equation}
    \frac{d\text{LeakyReLU}}{dz} = \begin{cases} 
    1 & \text{if } z > 0 \\
    \alpha & \text{if } z \leq 0
    \end{cases}
\end{equation}

where $\alpha$ is typically set to $0.01$. This small negative slope ensures gradient flow even for negative inputs.

\subsubsection{Linear Activation}

For the output layer in regression tasks, we use a linear (identity) activation:

\begin{equation}
    f(z) = z, \quad \frac{df}{dz} = 1
\end{equation}

This allows the network to produce unbounded continuous outputs suitable for regression problems.

\subsection{Cost Functions}

The cost function quantifies the discrepancy between predictions and true values, providing the optimization objective for training.

\subsubsection{Mean Squared Error (MSE)}

For regression tasks, we employ the MSE cost function:

\begin{equation}
    C(\mathbf{W}, \mathbf{b}) = \frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2
\end{equation}

where $n$ is the number of samples, $y_i$ are the true targets, and $\hat{y}_i$ are the predictions. The derivative with respect to predictions is:

\begin{equation}
    \frac{\partial C}{\partial \hat{y}_i} = \frac{2}{n}(\hat{y}_i - y_i)
\end{equation}

\subsubsection{Softmax Cross-Entropy}

For multiclass classification, we use the softmax cross-entropy loss. The softmax activation for the output layer computes class probabilities:

\begin{equation}
    \hat{y}_{ik} = \frac{e^{z_{ik}}}{\sum_{j=1}^{K} e^{z_{ij}}}
\end{equation}

where $K$ is the number of classes. The cross-entropy loss is:

\begin{equation}
    C(\mathbf{W}, \mathbf{b}) = -\frac{1}{n} \sum_{i=1}^{n} \sum_{k=1}^{K} y_{ik} \log(\hat{y}_{ik})
\end{equation}

where $y_{ik}$ is the one-hot encoded true label. The derivative with respect to logits has a remarkably simple form:

\begin{equation}
    \frac{\partial C}{\partial z_{ik}} = \frac{1}{n}(\hat{y}_{ik} - y_{ik})
\end{equation}

\subsection{Backpropagation Algorithm}

The backpropagation algorithm efficiently computes gradients of the cost function with respect to all network parameters using the chain rule of calculus. We present the algorithm for a general $L$-layer network.

\subsubsection{Forward Pass}

During forward propagation, we compute and store all intermediate activations:

\begin{align}
    \mathbf{a}^{(0)} &= \mathbf{X} \\
    \mathbf{z}^{(l)} &= \mathbf{W}^{(l)} \mathbf{a}^{(l-1)} + \mathbf{b}^{(l)}, \quad l = 1, \ldots, L \\
    \mathbf{a}^{(l)} &= f^{(l)}(\mathbf{z}^{(l)}), \quad l = 1, \ldots, L
\end{align}

\subsubsection{Backward Pass}

Starting from the output layer, we compute the error signal $\boldsymbol{\delta}^{(l)} = \frac{\partial C}{\partial \mathbf{z}^{(l)}}$ for each layer.

For the output layer $L$:
\begin{equation}
    \boldsymbol{\delta}^{(L)} = \frac{\partial C}{\partial \mathbf{a}^{(L)}} \odot f'^{(L)}(\mathbf{z}^{(L)})
\end{equation}

where $\odot$ denotes element-wise multiplication. For hidden layers $l = L-1, \ldots, 1$:

\begin{equation}
    \boldsymbol{\delta}^{(l)} = [(\mathbf{W}^{(l+1)})^T \boldsymbol{\delta}^{(l+1)}] \odot f'^{(l)}(\mathbf{z}^{(l)})
\end{equation}

The gradients with respect to weights and biases are then:

\begin{align}
    \frac{\partial C}{\partial \mathbf{W}^{(l)}} &= \frac{1}{n} (\mathbf{a}^{(l-1)})^T \boldsymbol{\delta}^{(l)} \\
    \frac{\partial C}{\partial \mathbf{b}^{(l)}} &= \frac{1}{n} \sum_{i=1}^{n} \boldsymbol{\delta}_i^{(l)}
\end{align}

where the averaging over $n$ samples accounts for batch processing.

\subsection{Optimization Methods}

We update network parameters using gradient-based optimization. The basic update rule for iteration $t$ is:

\begin{align}
    \mathbf{W}^{(l)}(t+1) &= \mathbf{W}^{(l)}(t) - \eta \frac{\partial C}{\partial \mathbf{W}^{(l)}} \\
    \mathbf{b}^{(l)}(t+1) &= \mathbf{b}^{(l)}(t) - \eta \frac{\partial C}{\partial \mathbf{b}^{(l)}}
\end{align}

where $\eta$ is the learning rate. We implemented and tested several optimization algorithms:

\subsubsection{Gradient Descent}

Standard gradient descent uses a fixed learning rate and computes gradients using the entire training set (or batches thereof). This serves as our baseline optimization method.

\subsubsection{RMSprop}

RMSprop adapts the learning rate for each parameter by maintaining a moving average of squared gradients:

\begin{align}
    v^{(l)}(t) &= \beta v^{(l)}(t-1) + (1-\beta) \left(\frac{\partial C}{\partial \mathbf{W}^{(l)}}\right)^2 \\
    \mathbf{W}^{(l)}(t+1) &= \mathbf{W}^{(l)}(t) - \frac{\eta}{\sqrt{v^{(l)}(t) + \epsilon}} \frac{\partial C}{\partial \mathbf{W}^{(l)}}
\end{align}

where $\beta \approx 0.9$ is the decay rate and $\epsilon = 10^{-8}$ prevents division by zero.

\subsubsection{Adam}

Adam (Adaptive Moment Estimation) combines momentum and adaptive learning rates by maintaining both first and second moment estimates:

\begin{align}
    m^{(l)}(t) &= \beta_1 m^{(l)}(t-1) + (1-\beta_1) \frac{\partial C}{\partial \mathbf{W}^{(l)}} \\
    v^{(l)}(t) &= \beta_2 v^{(l)}(t-1) + (1-\beta_2) \left(\frac{\partial C}{\partial \mathbf{W}^{(l)}}\right)^2 \\
    \hat{m}^{(l)}(t) &= \frac{m^{(l)}(t)}{1-\beta_1^t}, \quad \hat{v}^{(l)}(t) = \frac{v^{(l)}(t)}{1-\beta_2^t} \\
    \mathbf{W}^{(l)}(t+1) &= \mathbf{W}^{(l)}(t) - \frac{\eta}{\sqrt{\hat{v}^{(l)}(t)} + \epsilon} \hat{m}^{(l)}(t)
\end{align}

with typical values $\beta_1 = 0.9$ and $\beta_2 = 0.999$.

\subsection{Regularization}

To prevent overfitting, we can augment the cost function with $L_1$ or $L_2$ regularization penalties:

\begin{equation}
    C_{\text{reg}} = C + \lambda_2 \sum_{l} \|\mathbf{W}^{(l)}\|_F^2 + \lambda_1 \sum_{l} \|\mathbf{W}^{(l)}\|_1
\end{equation}

where $\|\mathbf{W}\|_F^2 = \sum_{ij} W_{ij}^2$ is the Frobenius norm and $\|\mathbf{W}\|_1 = \sum_{ij} |W_{ij}|$. These add corresponding terms to the gradients:

\begin{align}
    \frac{\partial C_{\text{reg}}}{\partial W_{ij}^{(l)}} &= \frac{\partial C}{\partial W_{ij}^{(l)}} + 2\lambda_2 W_{ij}^{(l)} + \lambda_1 \text{sign}(W_{ij}^{(l)})
\end{align}

\subsection{Data Preprocessing}

\subsubsection{Runge Function Dataset}

For regression tasks, we generated training data from the Runge function:

\begin{equation}
    f(x) = \frac{1}{1 + 25x^2}, \quad x \in [-1, 1]
\end{equation}

We sampled $n = 1000$ uniformly distributed points and optionally added Gaussian noise with standard deviation $\sigma = 0.1$. The data was split into training (70\%) and test (30\%) sets.

\subsubsection{MNIST Dataset}

For classification, we used the full MNIST dataset containing 70,000 grayscale images of handwritten digits (0-9), each of size $28 \times 28$ pixels. Each image was flattened into a 784-dimensional vector and normalized by dividing pixel values by 255 to obtain values in $[0, 1]$. The dataset was split into 60,000 training and 10,000 test samples following standard conventions.

\subsubsection{Feature Scaling}

For both tasks, we applied standardization to the input features:

\begin{equation}
    \tilde{x}_j = \frac{x_j - \mu_j}{\sigma_j}
\end{equation}

where $\mu_j$ and $\sigma_j$ are the mean and standard deviation of feature $j$ computed from the training set. This scaling improves convergence by ensuring all features have similar magnitudes.

\subsection{Implementation and Validation}

Our implementation follows an object-oriented design with modular components for activation functions, cost functions, and the neural network itself. The network architecture is specified by providing the input size, a list of hidden layer sizes, and corresponding activation functions.

To validate the correctness of our backpropagation implementation, we compared gradients against those computed using automatic differentiation via the \texttt{autograd} library. This serves as a numerical verification that our analytical gradient computation is correct.

We further validated our implementation by comparing performance against established machine learning libraries including Scikit-Learn, PyTorch, and TensorFlow/Keras for both regression and classification tasks.

\subsection{Performance Metrics}

For regression tasks, we use the coefficient of determination ($R^2$ score):

\begin{equation}
    R^2 = 1 - \frac{\sum_i (y_i - \hat{y}_i)^2}{\sum_i (y_i - \bar{y})^2}
\end{equation}

where $\bar{y}$ is the mean of true values. Additionally, we report MSE for direct comparison with Project 1 results.

For classification tasks, we use accuracy as the primary metric:

\begin{equation}
    \text{Accuracy} = \frac{1}{n}\sum_{i=1}^{n} \mathbb{I}(t_i = \arg\max_k \hat{y}_{ik})
\end{equation}

where $\mathbb{I}$ is the indicator function and $t_i$ is the true class label.

\subsection{Hyperparameter Exploration}

We systematically explored the hyperparameter space to identify optimal network architectures and training configurations. The explored parameters included:

\begin{itemize}
    \item Number of hidden layers: $\{0, 1, 2, 3\}$
    \item Nodes per hidden layer: $\{5, 10, 25, 50\}$
    \item Learning rates: $\eta \in \{0.001, 0.01, 0.1, 0.5, 1.0\}$
    \item Activation functions: Sigmoid, ReLU, Leaky ReLU
    \item Training epochs: $\{500, 1000, 1500, 2000\}$
    \item Regularization strengths: $\lambda \in \{0, 0.001, 0.01, 0.1\}$
\end{itemize}

Results are visualized using heatmaps to show the relationship between architectural choices and performance metrics, facilitating intuitive understanding of the design space.


\section{Results and discussion}
\subsection*{Regression}
To find the optimal architecture of our neural network for the regression task, we analyze the $R^2$-score of the model for different hyperparameters. The results are visualized as a heatmap in figure \ref{fig:nn_regression_heatmap}.
\begin{figure}[H]
    \centering
    \includegraphics[width=1\linewidth]{nn_regression_architecture_heatmap.pdf}
    \caption{Regression task. Neural network architecture comparison. 1500 epochs, 0.5 fixed learning rate gradient descent}
    \label{fig:nn_regression_heatmap}
\end{figure}
We can clearly see best performing architecture is the one with most layers and most nodes per layer.
\begin{figure}[H]
    \centering
    \includegraphics[width=1\linewidth]{nn_regression_predictions.pdf}
    \caption{Predictions of neural network on test data with optimal parameters found in heatmap: 3 hidden layers of 50 nodes each.}
    \label{fig:nn_regression_predictions}
\end{figure}
\textbf{Answer week 44 exercise 1c): }Figure \ref{fig:nn_regression_predictions} is unclear, as mathematical notation for the axis labels would be more precise. In addition to this the number of data points used and the fact that the prediction is performed on the test data (as opposed to the training data) should be clearly conveyed.

\subsection*{Classification}

\section{Conclusion}

% ===========================================
\onecolumngrid
% \bibliographystyle{apalike}
%\bibliographystyle{unsrt}
%\bibliography{ref}

\end{document}
