where $\prod_{j=k+1}^{t} \frac{\partial h_j}{\partial h_{j-1}}$ represents the chain of derivatives from time step $k$ to $t$. Unlike feedforward neural networks, where each input is processed independently, RNNs maintain a hidden state $h_t$ that acts as a dynamic memory, evolving recursively as the input sequence progresses. Formally, given an input sequence $\{x_1, x_2, \dots, x_T\}$, where $x_t \in \mathbb{R}^n$ represents the input vector at time $t$, the hidden state $h_t \in \mathbb{R}^m$ is updated via the recurrence relation:

$$h_t = f_h(W_{xh}x_t + W_{hh}h_{t-1} + b_h), \qquad (19.9)$$

where $W_{xh} \in \mathbb{R}^{m \times n}$, $W_{hh} \in \mathbb{R}^{m \times m}$, and $b_h \in \mathbb{R}^m$ are learnable parameters, and $f_h$ is a nonlinear activation function such as tanh or ReLU. The recursive structure inherently allows the hidden state $h_t$ to encode the entire history of the sequence up to time $t$. The output $y_t \in \mathbb{R}^o$ at each time step is computed as:

$$y_t = f_y(W_{hy}h_t + b_y), \qquad (19.10)$$

where $W_{hy} \in \mathbb{R}^{o \times m}$ and $b_y \in \mathbb{R}^o$ are additional learnable parameters, and $f_y$ is an optional output activation function, such as the softmax function for classification. To elucidate the recursive dynamics, we can expand $h_t$ explicitly in terms of the initial hidden state $h_0$ and all previous inputs $\{x_1, \dots, x_t\}$:

$$h_t = f_h(W_{xh}x_t + W_{hh}f_h(W_{xh}x_{t-1} + W_{hh}f_h(\dots f_h(W_{xh}x_1 + W_{hh}h_0 + b_h) + b_h) + b_h) + b_h). \quad (19.11)$$

This nested structure highlights the temporal dependencies and the potential challenges in training, such as the vanishing and exploding gradient problems. During training, the loss function $L$, which aggregates the discrepancies between the predicted outputs $y_t$ and the ground truth $y_t^{\text{true}}$, is typically defined as:

$$L = \sum_{t=1}^{T} \ell(y_t, y_t^{\text{true}}), \qquad (19.12)$$

where $\ell$ is a task-specific loss function, such as the mean squared error (MSE)

$$\ell(y, y^{\text{true}}) = \frac{1}{2} \|y - y^{\text{true}}\|^2 \qquad (19.13)$$

for regression or the cross-entropy loss for classification. To optimize $L$, gradient-based methods are employed, requiring the computation of derivatives of $L$ with respect to all parameters, such as $W_{xh}$, $W_{hh}$, and $b_h$. Using backpropagation through time (BPTT), the gradient of $L$ with respect to $W_{hh}$ is expressed as:

$$\frac{\partial L}{\partial W_{hh}} = \sum_{t=1}^{T} \sum_{k=1}^{t} \frac{\partial \ell_t}{\partial h_t} \prod_{j=k+1}^{t} \frac{\partial h_j}{\partial h_{j-1}} \frac{\partial h_k}{\partial W_{hh}}. \qquad (19.14)$$

Here,

$$\prod_{j=k+1}^{t} \frac{\partial h_j}{\partial h_{j-1}} \qquad (19.15)$$

is the product of Jacobian matrices that encode the influence of $h_k$ on $h_t$. The Jacobian $\frac{\partial h_j}{\partial h_{j-1}}$ itself is given by:

$$\frac{\partial h_j}{\partial h_{j-1}} = W_{hh} \odot f'_h(a_j), \qquad (19.16)$$

where

$$a_j = W_{xh}x_j + W_{hh}h_{j-1} + b_h, \qquad (19.17)$$

and $f'_h(a_j)$ denotes the elementwise derivative of the activation function. The repeated multiplication of these Jacobians can lead to exponential growth or decay of the gradients, depending on the spectral radius $\rho(W_{hh})$. Specifically, if $\rho(W_{hh}) > 1$, gradients explode, whereas if $\rho(W_{hh}) < 1$,
gradients vanish, severely hampering the training process for long sequences. To address these issues, modifications such as Long Short-Term Memory (LSTM) networks and Gated Recurrent Units (GRUs) introduce gating mechanisms that explicitly regulate the flow of information. In LSTMs, the cell state $c_t$, governed by additive dynamics, prevents vanishing gradients. The cell state is updated as:

$$c_t = f_t \odot c_{t-1} + i_t \odot \tanh(W_c x_t + U_c h_{t-1} + b_c), \quad (19.18)$$

where $f_t$ is the forget gate, $i_t$ is the input gate, and $U_c$, $W_c$, and $b_c$ are learnable parameters.

## 19.3 Sequence Modeling and Long Short-Term Memory (LSTM) and GRUs

### 19.3.1 Literature Review of Sequence Modeling and Long Short-Term Memory (LSTM) and GRUs

Potter and Egon (2024) [402] provided an extensive study of RNNs and their enhancements (LSTM and GRU) for time-series forecasting. The authors conduct an empirical comparison between these architectures and analyze their effectiveness in capturing long-term dependencies in sequential data. The study concludes that GRUs are computationally efficient but slightly less expressive than LSTMs, whereas standard RNNs suffer from vanishing gradients. Yatkin et. al. (2025) [403] introduced a topological perspective to RNNs, including LSTM and GRU, to address inconsistencies in real-world applications. The authors propose stability-enhancing mechanisms to improve RNN performance in finance and climate modeling. Their results show that topologically-optimized GRUs outperform traditional LSTMs in maintaining memory over long sequences. Saifullah (2024) [404] applied LSTM and GRU networks to biomedical image classification (chicken egg fertility detection). The paper demonstrates that GRU's simpler architecture leads to faster convergence while LSTMs achieve slightly higher accuracy due to better memory retention. The results highlight domain-specific strengths of LSTM vs. GRU, particularly in handling sparse feature representations. Alonso (2024) [405] rigorously explored the mathematical foundations of RNNs, LSTMs, and GRUs. The author provides a deep analysis of gating mechanisms, vanishing gradient solutions, and optimization techniques that improve sequence modeling. A theoretical comparison is drawn between hidden state dynamics in GRUs vs. LSTMs, supporting their application in NLP and time-series forecasting. Tu et. al. (2024) [406] in a medical AI study evaluates LSTMs and GRUs for predicting patient physiological metrics during sedation. The authors find that LSTMs retain more long-term dependencies in time-series medical data, making them suitable for patient monitoring, while GRUs are preferable for real-time predictions due to their lower computational overhead. Zuo et. al. (2025) [407] applied hybrid GRUs for predicting customer movements in stores using real-time location tracking. The authors propose a modified GRU-LSTM hybrid model that achieves state-of-the-art accuracy in trajectory prediction. The study demonstrates that GRUs alone may lack fine-grained memory retention, but a hybrid approach improves forecasting ability. Lima et. al. (2025) [408] developed an industrial AI application that demonstrated the efficiency of GRUs in process optimization. The study finds that GRUs outperform LSTMs in real-time predictive control of steel slab heating, showcasing their efficiency in applications where faster computations are required. Khan et. al. (2025) [409] integrated LSTMs with statistical ARIMA models to improve wind power forecasting. They demonstrate that hybrid LSTM-ARIMA models outperform standalone RNNs in handling weather-related sequential data, which is highly volatile. Guo and Feng (2024) [410] in an environmental AI study proposed a temporal attention-enhanced LSTM model to predict greenhouse climate variables. The research introduces a novel position-aware LSTM architecture that improves multi-step forecasting, which is critical for precision agriculture. Adegoke and Ghosh (2021) [1089] derived new infinite series that connect Fibonacci numbers with values of Riemann zeta functions via combinations of polygamma functions, thus illuminating deeper struc-
tural relations between these classical number-theoretic sequences. Guha and Ghosh (2021) [1110] introduce and systematically analyze the abundance index $In) = \frac{\sigma(n)}{n}$, explore properties of its value set, and propose a novel class of superabundant numbers, investigating their connections to the Riemann Hypothesis. Adegoke (2021) et. al. [1117] evaluated some cubic binomial Fibonacci and Lucas sums. Bruun & Ghosh (2021) [1149] define an adjacency-matrix representation of the Collatz graph and show how a flow-diagram “picture” arises from it, arguing that the $3n+1$ problem can be viewed in terms of groups of values rather than isolated trajectories, and indicating that knowing the stopping time for one element of a group may allow prediction for the entire group. Ghosh (2020) [1163] provides a comprehensive overview of lattices and their applications in number theory, particularly in the study of Diophantine equations and functional analysis, highlighting their significance in the geometry of numbers. Cox and Ghosh (2022) [1191] introduce an analogue of Lagarias’ inequality for the Riemann Hypothesis, utilizing a function with properties similar to those of harmonic numbers, thereby extending the scope of inequalities pertinent to the hypothesis. Cox et. al. (2022) [1202] investigate the relationship between abundant numbers and the Riemann Hypothesis, providing insights into how certain classes of numbers may influence the hypothesis. Cox et. al. (2021) [1211] explore the implications of Fermat’s Last Theorem and related problems, providing new insights into the connections between number theory and algebraic structures. Abdelhamid (2024) [411] explored IoT-based energy forecasting using deep RNN architectures, including LSTM and GRU. The study concludes that GRUs provide faster inference speeds but LSTMs capture more accurate long-range dependencies, making them more reliable for complex forecasting.

## 19.3.2 Analysis of Sequence Modeling and Long Short-Term Memory (LSTM) and GRUs

Sequence modeling in Recurrent Neural Networks (RNNs) represents a powerful framework for capturing temporal dependencies in sequential data, enabling the learning of both short-term and long-term patterns. The primary characteristic of RNNs lies in their recurrent architecture, where the hidden state $h_t$ at time step $t$ is updated as a function of both the current input $x_t$ and the hidden state at the previous time step $h_{t-1}$. Mathematically, this recurrent relationship can be expressed as:

$$h_t = f(W_h h_{t-1} + W_x x_t + b_h) \qquad (19.19)$$

Here, $W_h$ and $W_x$ are weight matrices corresponding to the previous hidden state $h_{t-1}$ and the current input $x_t$, respectively, while $b_h$ is a bias term. The function $f(\cdot)$ is a non-linear activation function, typically chosen as the hyperbolic tangent tanh or rectified linear unit (ReLU). The output $y_t$ at each time step is derived from the hidden state $h_t$ through a linear transformation, followed by a non-linear activation, yielding:

$$y_t = g(W_y h_t + b_y) \qquad (19.20)$$

where $W_y$ is the weight matrix connecting the hidden state to the output space, and $b_y$ is the associated bias term. The function $g(\cdot)$ is generally a softmax activation for classification tasks or a linear activation for regression problems. The key feature of this structure is the interdependence of the hidden state across time steps, allowing the model to capture the history of past inputs and produce predictions that incorporate temporal context. Training an RNN involves minimizing a loss function $L$, which quantifies the discrepancy between the predicted outputs $y_t$ and the true outputs $y_t^{\text{true}}$ across all time steps. A common loss function used in classification tasks is the cross-entropy loss, while regression tasks often utilize mean squared error. To optimize the parameters of the network, the model employs **Backpropagation Through Time (BPTT)**, a variant of the standard backpropagation algorithm adapted for sequential data. The primary challenge in BPTT arises from the recurrent nature of the network, where the hidden state at each time step depends on the previous hidden state. The gradient of the loss function with respect to the hidden state at
time step $t$ is computed recursively, reflecting the temporal structure of the model. The chain rule is applied to compute the gradient of the loss with respect to the hidden state:

$$ \frac{\partial L}{\partial h_t} = \frac{\partial L}{\partial y_t} \cdot \frac{\partial y_t}{\partial h_t} + \sum_{t'=t+1}^{T} \frac{\partial L}{\partial h_{t'}} \cdot \frac{\partial h_{t'}}{\partial h_t} \qquad (19.21) $$

Here, $\frac{\partial L}{\partial y_t}$ is the gradient of the loss with respect to the output, and $\frac{\partial y_t}{\partial h_t}$ represents the Jacobian of the output with respect to the hidden state. The second term in this expression corresponds to the accumulated gradients propagated from future time steps, incorporating the temporal dependencies across the entire sequence. This recursive gradient calculation allows for updating the weights and biases of the RNN, adjusting them to minimize the total error across the sequence. The gradients of the loss function with respect to the parameters of the network, such as $W_h, W_x$, and $W_y$, are computed using the chain rule. For example, the gradient of the loss with respect to $W_x$ is:

$$ \frac{\partial L}{\partial W_x} = \sum_{t=1}^{T} \frac{\partial L}{\partial h_t} \cdot x_t^\top \qquad (19.22) $$

This captures the contribution of each input to the overall error at all time steps, ensuring that the model learns the correct relationships between inputs and hidden states. Similarly, the gradients with respect to $W_h$ and $b_h$ account for the recurrence in the hidden state, enabling the model to adjust its internal parameters in response to the sequential nature of the data. Despite their theoretical elegance, RNNs face significant practical challenges during training, primarily due to the **vanishing gradients problem**. This issue arises when the gradients propagate through many time steps, causing them to decay exponentially, especially when using activation functions like tanh. As a result, the influence of distant time steps diminishes, making it difficult for the network to learn long-term dependencies. The mathematical manifestation of this problem is seen in the norm of the Jacobian matrices associated with the hidden state updates. If the spectral radius of the weight matrices $W_h$ is close to or greater than 1, the gradients can either vanish or explode, leading to unstable training dynamics. To mitigate this issue, several solutions have been proposed, including the use of Long Short-Term Memory (LSTM) networks and Gated Recurrent Units (GRUs), which introduce gating mechanisms to better control the flow of information through the network. LSTMs, for example, incorporate a memory cell $C_t$, which allows the network to store information over long periods of time. The update rules for the LSTM are governed by three gates: the forget gate $f_t$, the input gate $i_t$, and the output gate $o_t$, which control how much of the previous memory and new information to retain. The equations governing the LSTM are:

$$ f_t = \sigma(W_f h_{t-1} + U_f x_t + b_f) \qquad (19.23) $$

$$ i_t = \sigma(W_i h_{t-1} + U_i x_t + b_i) \qquad (19.24) $$

$$ o_t = \sigma(W_o h_{t-1} + U_o x_t + b_o) \qquad (19.25) $$

$$ \tilde{C}_t = \tanh(W_C h_{t-1} + U_C x_t + b_C) \qquad (19.26) $$

$$ C_t = f_t \cdot C_{t-1} + i_t \cdot \tilde{C}_t \qquad (19.27) $$

$$ h_t = o_t \cdot \tanh(C_t) \qquad (19.28) $$

In these equations, the forget gate $f_t$ determines how much of the previous memory cell $C_{t-1}$ to retain, the input gate $i_t$ governs how much new information to store in the candidate memory cell $\tilde{C}_t$, and the output gate $o_t$ controls how much of the memory cell should influence the current output. The LSTM's architecture allows for the maintenance of long-term dependencies by selectively forgetting or retaining information, effectively alleviating the vanishing gradient problem and enabling the network to learn from longer sequences. The GRU, an alternative to the LSTM, simplifies this architecture by combining the forget and input gates into a single update gate $z_t$,
and introduces a reset gate $r_t$ to control the influence of the previous hidden state. The GRU's update rules are:

$$z_t = \sigma(W_z h_{t-1} + U_z x_t + b_z) \qquad (19.29)$$

$$r_t = \sigma(W_r h_{t-1} + U_r x_t + b_r) \qquad (19.30)$$

$$\tilde{h}_t = \tanh(Wh_{t-1} + Ux_t + b) \qquad (19.31)$$

$$h_t = (1 - z_t) \cdot h_{t-1} + z_t \cdot \tilde{h}_t \qquad (19.32)$$

Here, $z_t$ controls the amount of the previous hidden state $h_{t-1}$ to retain, and $r_t$ determines how much of the previous hidden state should influence the candidate hidden state $\tilde{h}_t$. The GRU's simplified structure still allows it to effectively capture long-range dependencies while being computationally more efficient than the LSTM.

In summary, sequence modeling in RNNs involves a series of recurrent updates to the hidden state, driven by both the current input and the previous hidden state, and is trained via backpropagation through time. The introduction of specialized gating mechanisms in LSTMs and GRUs alleviates issues such as vanishing gradients, enabling the networks to learn and maintain long-term dependencies. Through these advanced architectures, RNNs can effectively model complex temporal relationships, making them powerful tools for tasks such as time-series prediction, natural language processing, and sequence generation.

# 19.4 Applications in Natural Language Processing

## 19.4.1 Literature Review of Applications in Natural Language Processing

Yang et. al. (2020) [392] explored the effectiveness of deep learning models, including RNNs, for sentiment analysis in e-commerce platforms. It emphasizes how RNN architectures, including LSTMs and GRUs, outperform traditional NLP techniques by capturing sequential dependencies in customer reviews. The study provides empirical evidence demonstrating the superior accuracy of RNNs in analyzing consumer sentiment. Manikandan et. al. (2025) [393] investigated how RNNs can improve spam detection in email filtering. By leveraging recurrent structures, the study demonstrates how RNNs effectively identify patterns in email text that indicate spam or phishing attempts. It also compares RNN-based models with other ML approaches, highlighting the robustness of RNNs in handling contextual word sequences. Isiaka et. al. (2025) [394] examined AI technologies, particularly deep learning models, for predictive healthcare applications. It highlights how RNNs can analyze patient records and medical reports using NLP techniques. The study shows that RNN-based NLP models enhance medical diagnostics and decision-making by extracting meaningful insights from unstructured text data. Petrov et. al. (2025) [395] discussed the role of RNNs in emotion classification from textual data, an essential NLP task. The paper evaluates various RNN-based architectures, including BiLSTMs, to enhance the accuracy of emotion recognition in social media texts and chatbot responses. Liang (2025) [396] focused on the application of RNNs in educational settings, specifically for automated grading and feedback generation. The study presents an RNN-based NLP system capable of analyzing student responses, providing real-time assessments, and generating contextual feedback. Jin (2025) [397] explored how RNNs optimize text generation tasks related to pharmaceutical education. It demonstrates how NLP-powered RNN models generate high-quality textual summaries from medical literature, ensuring accurate knowledge dissemination in the pharmaceutical industry. McNicholas et. al. (2025) [398] investigated how RNNs facilitate clinical decision-making in critical care by extracting insights from unstructured medical text. The research highlights how RNN-based NLP models enhance patient care by predicting outcomes based on clinical notes and physician reports. Abbas
and Khammas (2024) [399] introduced an RNN-based NLP technique for detecting malware in IoT networks. The study illustrates how RNN classifiers process logs and textual patterns to identify malicious software, making RNNs crucial in cybersecurity applications. Kalonia and Upadhyay (2025) [400] applied RNNs to software fault prediction using NLP techniques. It shows how recurrent networks analyze bug reports and software documentation to predict potential failures in software applications, aiding developers in proactive debugging. Cox et. al. (2021) [1220] explore the relationship between the Farey sequence and the Mertens function, introducing functions analogous to the Mertens function for subsets of fractions in Farey sequences and proposing a conjecture linking the Farey sequence to the Riemann Hypothesis. Jormakka and Ghosh [1229] provide a new proof demonstrating that a prime $d \equiv 3 \pmod 8$ is not a congruent number, utilizing a method that can be generalized to determine the congruency of non-prime numbers without relying on the Birch and Swinnerton-Dyer conjecture. Cox et. al. (2021) [1239] investigate the sign behavior of the Mertens function by introducing functions analogous to the Mertens function for subsets of fractions in Farey sequences, establishing upper bounds for functions involving the sum of the signs of Mertens function values, and exploring relationships between the Farey sequence and the Riemann hypothesis. Cox et. al. (2021) [1247] analyze the energy spectral density of the Mertens function, demonstrating that the sum of squared values of the Mertens function over Farey sequences provides upper bounds that support the Riemann Hypothesis. Cox and Ghosh (2021) [1255] extended the Prime Number Theorem by introducing a generalized asymptotic formula for the sum of divisors function, providing deeper insights into the distribution of divisors and their connection to prime numbers. Cox et. al. (2020) [1261] establish a novel connection between nth-degree polynomial congruences modulo a prime (with n roots) and higher-order Fibonacci and Lucas recurrences, deriving new relations between the roots of recurrence sequences and polynomial congruence behavior. Maji et. al. (2021) [1271] test the conjecture of equivalence of nonequilibrium ensembles in two-dimensional turbulence by constructing a formally time-reversible Navier-Stokes system with conserved energy and enstrophy, and show that its statistical properties closely match those of the standard (irreversible) Navier-Stokes system across both inverse and direct cascade regimes. Ghosh (2022) [1277] offers a new proof of the irrationality of $N^{p/q}$ by leveraging the binomial theorem in a novel way, enriching the collection of elementary methods in transcendence theory. Han et. al. (2025) [401] discussed RNN applications in conversational AI, focusing on chatbots and virtual assistants. The study presents an RNN-driven NLP model for improving dialogue management and user interactions, significantly enhancing the responsiveness of AI-powered chat systems.

## 19.4.2 Analysis of Applications in Natural Language Processing

Recurrent Neural Networks (RNNs) are deep learning architectures that are explicitly designed to handle sequential data, a key feature that makes them indispensable for applications in Natural Language Processing (NLP). The mathematical foundation of RNNs lies in their ability to process sequences of inputs, $x_1, x_2, \dots, x_T$, where T denotes the length of the sequence. At each time step t, the network updates its hidden state, $h_t$, using both the current input $x_t$ and the previous hidden state $h_{t-1}$. This recursive relationship is represented mathematically by the following equation:

$$h_t = \sigma(W_h h_{t-1} + W_x x_t + b) \qquad (19.33)$$

Here, $\sigma$ is a nonlinear activation function such as the hyperbolic tangent (tanh) or the rectified linear unit (ReLU), $W_h$ is the weight matrix associated with the previous hidden state $h_{t-1}$, $W_x$ is the weight matrix associated with the current input $x_t$, and $b$ is a bias term. The nonlinearity introduced by $\sigma$ allows the network to learn complex relationships between the input and the output. The output $y_t$ at each time step is computed as:

$$y_t = W_y h_t + c \qquad (19.34)$$

where $W_y$ is the weight matrix corresponding to the output and $c$ is the bias term for the output. The output $y_t$ is then used to compute the predicted probability distribution over possible outputs
at each time step, typically through a softmax function for classification tasks:

$$P(y_t|h_t) = \mathrm{softmax}(W_y h_t + c) \qquad (19.35)$$

In NLP tasks such as language modeling, the objective is to predict the next word in a sequence given all previous words. The RNN is trained to estimate the conditional probability distribution $P(w_t|w_1, w_2, \dots, w_{t-1})$ of the next word $w_t$ based on the previous words. The full likelihood of the sequence $w_1, w_2, \dots, w_T$ can be written as:

$$P(w_1, w_2, \dots, w_T) = \prod_{t=1}^{T} P(w_t|w_1, w_2, \dots, w_{t-1}) \qquad (19.36)$$

For an RNN, this conditional probability is modeled by recursively updating the hidden state and generating a probability distribution for each word. At each time step, the probability of the next word is computed as:

$$P(w_t|h_{t-1}) = \mathrm{softmax}(W_y h_t + c) \qquad (19.37)$$

The network is trained by minimizing the negative log-likelihood of the true word sequence:

$$\mathcal{L} = - \sum_{t=1}^{T} \log P(w_t|h_{t-1}) \qquad (19.38)$$

This loss function guides the optimization of the weight matrices $W_h$, $W_x$, and $W_y$ to maximize the likelihood of the correct word sequences. As the network learns from large datasets, it develops the ability to predict words based on the context provided by previous words in the sequence. A key extension of RNNs in NLP is machine translation, where one sequence of words in one language is mapped to another sequence in a target language. This is typically modeled using sequence-to-sequence (Seq2Seq) architectures, which consist of two RNNs: the encoder and the decoder. The encoder RNN processes the input sequence $x_1, x_2, \dots, x_T$, updating its hidden state at each time step:

$$h_t^{\mathrm{enc}} = \sigma(W_h^{\mathrm{enc}} h_{t-1}^{\mathrm{enc}} + W_x^{\mathrm{enc}} x_t + b^{\mathrm{enc}}) \qquad (19.39)$$

The final hidden state $h_T^{\mathrm{enc}}$ of the encoder is passed to the decoder as its initial hidden state. The decoder RNN generates the target sequence $y_1, y_2, \dots, y_T$ by updating its hidden state at each time step, using both the previous hidden state $h_{t-1}^{\mathrm{dec}}$ and the previous output $y_{t-1}$:

$$h_t^{\mathrm{dec}} = \sigma(W_h^{\mathrm{dec}} h_{t-1}^{\mathrm{dec}} + W_x^{\mathrm{dec}} y_{t-1} + b^{\mathrm{dec}}) \qquad (19.40)$$

The decoder produces a probability distribution over the target vocabulary at each time step:

$$P(y_t|h_t^{\mathrm{dec}}) = \mathrm{softmax}(W_y^{\mathrm{dec}} h_t^{\mathrm{dec}} + c^{\mathrm{dec}}) \qquad (19.41)$$

The training of the Seq2Seq model is based on minimizing the cross-entropy loss function:

$$\mathcal{L} = - \sum_{t=1}^{T} \log P(y_t|h_t^{\mathrm{dec}}) \qquad (19.42)$$

This ensures that the network learns to map input sequences to output sequences. By training on a large corpus of paired sentences, the Seq2Seq model learns to translate sentences from one language to another, with the encoder capturing the context of the input sentence and the decoder generating the translated sentence.

RNNs are also effective in sentiment analysis, a task where the goal is to classify the sentiment
of a sentence (positive, negative, or neutral). Given a sequence of words $x_1, x_2, \dots, x_T$, the RNN processes each word sequentially, updating its hidden state:

$$h_t = \sigma(W_h h_{t-1} + W_x x_t + b) \qquad (19.43)$$

After processing the entire sentence, the final hidden state $h_T$ is used to classify the sentiment. The output is obtained by applying a softmax function to the final hidden state:

$$y = \mathrm{softmax}(W_y h_T + c) \qquad (19.44)$$

where $W_y$ is the weight matrix associated with the output layer. The network is trained to minimize the cross-entropy loss:

$$\mathcal{L} = - \sum_{t=1}^{T} \log P(y|h_T) \qquad (19.45)$$

This allows the RNN to classify the overall sentiment of the sentence by learning the relationships between words and sentiment labels. Sentiment analysis is useful for applications such as customer feedback analysis, social media monitoring, and opinion mining. In Named Entity Recognition (NER), RNNs are used to identify and classify named entities, such as people, locations, and organizations, in a text. The RNN processes each word $x_t$ in the sequence, updating its hidden state at each time step:

$$h_t = \sigma(W_h h_{t-1} + W_x x_t + b) \qquad (19.46)$$

The output at each time step is a probability distribution over possible entity labels:

$$P(y_t|h_t) = \mathrm{softmax}(W_y h_t + c) \qquad (19.47)$$

The network is trained to minimize the cross-entropy loss:

$$\mathcal{L} = - \sum_{t=1}^{T} \log P(y_t|h_t) \qquad (19.48)$$

By learning to classify each word with the appropriate entity label, the RNN can perform information extraction tasks, such as identifying people, organizations, and locations in text. This is crucial for applications such as document categorization, knowledge graph construction, and question answering. In speech recognition, RNNs are used to transcribe spoken language into text. The input to the RNN consists of a sequence of acoustic features, such as Mel-frequency cepstral coefficients (MFCCs), which are extracted from the audio signal. At each time step $t$, the RNN updates its hidden state:

$$h_t = \sigma(W_h h_{t-1} + W_x x_t + b) \qquad (19.49)$$

The output at each time step is a probability distribution over phonemes or words:

$$P(w_t|h_t) = \mathrm{softmax}(W_y h_t + c) \qquad (19.50)$$

The network is trained by minimizing the negative log-likelihood:

$$\mathcal{L} = - \sum_{t=1}^{T} \log P(w_t|h_t) \qquad (19.51)$$

By learning the mapping between acoustic features and corresponding words or phonemes, the RNN can transcribe speech into text, which is fundamental for applications such as voice assistants, transcription services, and speech-to-text systems.

In summary, RNNs are powerful tools for processing sequential data in NLP tasks such as machine translation, sentiment analysis, named entity recognition, and speech recognition. Their ability to process input sequences in a time-dependent manner allows them to capture long-range dependencies, making them well-suited for complex tasks in NLP and beyond. However, challenges such as the vanishing and exploding gradient problems necessitate the use of more advanced architectures, like LSTMs and GRUs, to enhance their performance in real-world applications.
## 19.5 Deep Learning and the Collatz Conjecture

The **Collatz conjecture**, also known as the *3x + 1 problem*, is one of the most famous unsolved problems in number theory. It is defined by the recurrence relation:

$$ T(n) = \begin{cases} \frac{n}{2}, & \text{if } n \text{ is even} \\ 3n + 1, & \text{if } n \text{ is odd} \end{cases} \tag{19.52} $$

for any positive integer $n$. The conjecture asserts that for all $n \in \mathbb{N}$, repeated iteration of $T(n)$ eventually results in 1, after which the sequence enters the cycle $1 \to 4 \to 2 \to 1$. Despite its simple formulation, proving this statement for all $n$ has remained elusive.

### 19.5.1 Literature Review of Deep Learning and the Collatz Conjecture

Several approaches attempt to model the behavior of the Collatz sequence using transformations and modular arithmetic. Notably, Rahn et. al. (2021) [915] proposed an algorithm that linearizes Collatz convergence by introducing arithmetic rules that govern the entire dynamical system. This approach provides a structured way to analyze sequence behavior, which can be beneficial in training deep learning models for discrete dynamical problems. Cox et. al. (2021) [722] extended the analysis of the classic Collatz conjecture to generalized sequences defined by the recurrence relation $3n + c$, where $c$ is an odd integer not divisible by 3. The authors build upon techniques introduced by Halbeisen and Hungerbühler to analyze rational Collatz cycles, applying these methods to the more general $3n + c$ sequences. In these generalized sequences, cycles are categorized based on their length and the number of odd elements they contain. The study introduces the concept of parity vectors—sequences of 0s and 1s representing even and odd elements, respectively—to determine the smallest and largest odd elements within these cycles. The authors observe that these smallest and largest odd elements often belong to the same cycle, suggesting that the parity vector generated by the floor function can be rotated to match that generated by the ceiling function. This rotation involves two linear congruences, and the natural numbers produced by one of these congruences appear to be uniformly distributed after sorting, exhibiting properties similar to the zeros of the Riemann zeta function. The paper provides a detailed mathematical framework for analyzing these cycles, offering insights into the behavior of generalized $3n + c$ sequences and their relation to classical problems in number theory. Jormakka and Ghosh (2021) [1286] model the persistence of recessive lethal alleles in populations and show that such alleles can remain over long time scales if certain family lineages have disproportionately large numbers of offspring, thus providing an alternative explanation to heterozygote advantage or genetic drift. Jormakka and Ghosh (2021) [1297] develop a Markov Decision Process model for intruder detection scenarios, analyze the associated risks and costs under various strategies, and derive optimal policies balancing detection efficacy with operational expense. Jormakka and Ghosh (2021) [1298] derive methods to approximate and more exactly solve the Howard equation in multiservice loss systems, enabling them to compute not only the expected average cost under congestion pricing but also the full distribution of costs (i.e. the risk that actual costs deviate from the mean). Krishna et. al. (2021) [1299] introduce and analyze the concept of *special primes*, characterizing them by the condition

$$ q \equiv 1 \pmod{p} \quad \text{and} \quad 2^{\frac{q-1}{p}} \equiv 1 \pmod{q}, $$

proving that $q$ is a special prime of $p$ exactly when the order of 2 modulo $q$ divides $\frac{q-1}{p}$, and demonstrating existence results and negative criteria based on generator properties of the multiplicative group. Gupta et. al. (2022) [1312] analyze *rotating binaries* by proving that in certain cases a left-rotation of a binary number by $r$ bits yields the divisibility

$$ l \mid N_1(B) \cdot r + 1, $$
where $l$ is its bit length and $N_1(B)$ its Hamming weight, and they extend this result to generalized $kn+c$ cycles, narrowing the conditions under which such divisibility holds. Cox et. al. (2022) [1323] introduce a unified framework for analyzing reciprocity laws across multiple degrees, from quadratic to quintic, by employing group-theoretic methods to determine the 'degrees of freedom' in residue classes modulo primes, thereby extending classical reciprocity to higher powers and providing new insights into the structure of finite abelian groups. Cox and Ghosh (2023) [1338] introduce a variant of the Mertens function with exclusively negative values (except at $M(1)$), explore its relationship with the sum of divisors function, and investigate its potential applications in partially factoring certain natural numbers. Ghosh (2023) [1433] demonstrates the equivalence between the Galerkin finite element and central difference formulations of the advection-diffusion equation by deriving a unified discretization scheme that preserves both mass and momentum conservation, and provides a framework for comparing the stability and accuracy of both methods. Ghosh (2023) [1434] conducts a comprehensive stability analysis of three time integration schemes—Forward Euler, Backward Euler, and Central Difference—by deriving their respective stability regions and providing detailed comparisons to guide their application in solving partial differential equations.

### 19.5.2 Analysis of Deep Learning and the Collatz Conjecture

Since the transformation is entirely deterministic, a neural network would not necessarily need to approximate probabilities but rather learn the structural properties of the sequence. Deep learning techniques such as recurrent neural networks (RNNs), transformers, and graph neural networks (GNNs) provide a potential approach to uncovering hidden patterns that could offer insights into this problem. From a machine learning perspective, the Collatz sequence can be framed as a sequence prediction problem, where the goal is to predict the next number in the sequence given the previous values. Given an initial number $n_0$, the sequence

$$ (n_0, n_1, n_2, \dots, n_k) \qquad (19.53) $$

can be viewed as a time-series dataset. A Recurrent Neural Network (RNN), such as a Long Short-Term Memory (LSTM) network or a Transformer model, can be trained on known Collatz sequences to generate sequences for new unseen numbers. The function

$$ \hat{T}_{\theta}(n) = f_{\theta}(n) \qquad (19.54) $$

where $f_{\theta}$ is a neural network parameterized by $\theta$, can be trained to approximate $T(n)$. The loss function for training would be:

$$ \mathcal{L}(\theta) = \sum_{n=1}^{N} |f_{\theta}(n) - T(n)|^2. \qquad (19.55) $$

Another perspective on the Collatz problem is to consider it as a **graph learning task**. The transformation $T(n)$ induces a **directed graph** $G(V, E)$, where each node $V$ represents an integer $n$ and each directed edge represents the transformation $n \to T(n)$. A **Graph Neural Network (GNN)** can be trained on this graph structure to embed numerical relationships into a **low-dimensional vector space**, revealing hidden clustering patterns. The message-passing function in a GNN is given by:

$$ h_{n}^{(t+1)} = \sigma \left( W \sum_{m \in \mathcal{N}(n)} h_{m}^{(t)} \right), \qquad (19.56) $$

where $h_n^{(t)}$ is the embedding of node $n$ at iteration $t$, $W$ is a trainable weight matrix, and $\mathcal{N}(n)$ represents the set of neighboring nodes. Even though the Collatz function is deterministic, we can
model it as a **stochastic process** by introducing a probabilistic relaxation. Consider a Markov process $X_t$ defined on $\mathbb{N}$ where the transition probabilities are:

$$P(X_{t+1} = m | X_t = n) = \alpha\delta(m - T(n)) + (1 - \alpha)p(m), \qquad (19.57)$$

where $\delta(\cdot)$ is the Dirac delta function and $p(m)$ is a small noise term. A **reinforcement learning (RL)** agent can be defined with a policy $\pi(n)$ that selects transitions based on a learned strategy. The goal is to minimize the cost function:

$$J(\pi) = \mathbb{E}_{\pi} \left[ \sum_{t=0}^{\infty} \gamma^t C(n_t) \right], \qquad (19.58)$$

where $C(n_t)$ is a cost function that represents the number of steps needed to reach 1. RL-based approaches can be used to estimate statistical properties of the stopping time function $S(n)$, potentially revealing asymptotic behaviors in large $n$. Deep learning can serve as a heuristic tool to analyze large-scale computational data, offering numerical insights that might guide future theoretical approaches. If certain statistical properties—such as mean stopping times or clustering of trajectory lengths—can be detected via deep learning, they could provide new conjectures or constraints that mathematicians can then attempt to prove rigorously.

The intersection of deep learning and the Collatz conjecture is an exciting and emerging field. While deep learning cannot replace traditional mathematical proofs, it provides a powerful experimental tool to analyze the statistical properties of the Collatz sequence, detect hidden structures in its graph representation, and study the problem through reinforcement learning and Markov models. By leveraging graph neural networks, recurrent models, and probabilistic relaxations, we can explore new perspectives on this classic problem. However, the deterministic nature of the Collatz function suggests that deep learning alone will not yield a proof, but rather assist in formulating new hypotheses and guiding analytical techniques.

## 19.6 Mertens Function and the Collatz Conjecture

The **Mertens function** is defined as

$$M(n) = \sum_{k=1}^{n} \mu(k), \qquad (19.59)$$

where the Möbius function $\mu(n)$ satisfies

$$\mu(n) = \begin{cases} 1, & \text{if } n \text{ is square-free with an even number of prime factors,} \\ -1, & \text{if } n \text{ is square-free with an odd number of prime factors,} \\ 0, & \text{if } n \text{ is divisible by a square } p^2 \text{ for some prime } p. \end{cases} \qquad (19.60)$$

The asymptotic growth of $M(n)$ has a deep connection to the **Riemann Hypothesis (RH)**. The classical Mertens conjecture proposed

$$|M(n)| \le \sqrt{n} \quad \forall n \ge 1, \qquad (19.61)$$

which, if true, would have implied RH. However, it was disproven by Odlyzko and te Riele (1985). The best known upper bound under RH is

$$M(n) = O(n^{\frac{1}{2}+\epsilon}) \quad \text{for any } \epsilon > 0. \qquad (19.62)$$
Without assuming RH, the best unconditional bound known is

$$M(n) = O\left(\frac{n}{\exp(c\sqrt{\log n})}\right), \quad c > 0. \qquad (19.63)$$

A fundamental connection between $M(n)$ and the **zeros of the Riemann zeta function** arises via

$$\zeta(s) = \sum_{n=1}^{\infty} \frac{1}{n^s}, \quad \text{Re}(s) > 1. \qquad (19.64)$$

By analytic continuation, $\zeta(s)$ satisfies the **functional equation**

$$\zeta(s) = 2^s \pi^{s-1} \sin\left(\frac{\pi s}{2}\right) \Gamma(1-s)\zeta(1-s). \qquad (19.65)$$

## 19.6.1 Literature Review of Bounds of Mertens Function

Let us discuss the literature review on bounds for the Mertens function $M(x) = \sum_{n \le x} \mu(n)$ (where $\mu$ is the Möbius function). We summarize major historical milestones, unconditional vs. conditional (Riemann-hypothesis) results, computational advances that give explicit numeric bounds, the main techniques used, and current open problems — with pointers to the primary sources for each load-bearing claim.

The object of study.

$$M(x) := \sum_{n \le x} \mu(n) \qquad (19.66)$$

and the classical Mertens conjecture asserted $|M(x)| \le \sqrt{x}$ for all $x$. This conjecture tied into the distribution of zeta zeros because smallness of $M(x)$ is a strong sieve-type statement about primes and zeros.

**Historical turning point — disproof of the Mertens conjecture:** Odlyzko and the Riele (1985) [865] used lattice-reduction ideas (LLL and its use in a diophantine approximation translation of the “large oscillation” problem) to prove the Mertens conjecture false: they showed

$$\limsup_{x \to \infty} \frac{M(x)}{\sqrt{x}} > 1.06, \quad \liminf_{x \to \infty} \frac{M(x)}{\sqrt{x}} < -1.009, \qquad (19.67)$$

so $|M(x)| > \sqrt{x}$ occurs infinitely often (existence of a counterexample is guaranteed), although an explicit small $x$ counterexample was not found. Their paper explains the reduction to (simultaneous weighted) inhomogeneous diophantine approximation and the role of LLL.

**Unconditional analytic bounds (classical results):** Even long before and apart from Mertens' conjecture, analytic number theory produced unconditional upper bounds that are far weaker than $\sqrt{x}$ but still show $M(x)$ is small compared with $x$. A standard classical result by Davenport (1937) [866] is that for every fixed $A > 0$,

$$M(x) = O\left(\frac{x}{(\log x)^A}\right), \qquad (19.68)$$

so $M(x) = o(x)$ and stronger logarithmic savings hold by classical multiplicative-function methods. These unconditional bounds come from mean-value and exponential-sum estimates and follow from the prime number theorem + Fourier methods.

**Conditional (RH/GRH) bounds and modern refinements:** The Riemann Hypothesis (RH)
controls the size of $M(x)$ very tightly. Equivalently one may state (standard fact) that RH is equivalent to

$$M(x) = O(x^{1/2+\epsilon}) \quad \text{for every } \epsilon > 0, \qquad (19.69)$$

and more refined conditional upper bounds have been proved assuming RH (and sometimes additional hypotheses about zeta moments or zero spacing). The state-of-the-art conditional result is due to Soundararajan (2007) [867] & Soundararajan (2009) [868]: assuming RH he obtains an upper bound of the form

$$M(x) \ll \sqrt{x}, \exp!((\log x)^{1/2}(\log \log x)^C) \qquad (19.70)$$

for an explicit constant $C$; the method gives the best known improvement on older Titchmarsh/Hardy-Littlewood type bounds. Maier-Montgomery and other authors developed alternative conditional bounds based on zero-spacing and large deviations of zero counts; Ng (2004) [1076] and others studied distributional consequences (limiting distributions of normalized $M(e^y)$) under RH + additional moment conjectures.

**Probabilistic / heuristic conjectures about the true order:** Heuristics from random-matrix / probabilistic modeling of zeta zeros suggest that the true maximal order of $(M(x))$ is slightly larger than $\sqrt{x}$ by iterated logarithmic factors. For example, (unpublished) Gonek conjectured that

$$\limsup_{x \to \infty} \frac{|M(x)|}{\sqrt{x}(\log \log \log x)^{5/4}} \in (0, \infty), \qquad (19.71)$$

and Nathan Ng gives probabilistic evidence (conditional on RH + conjectures on zeta negative moments) supporting such triple-log factors and shows the existence of a limiting distribution for appropriately scaled $M(e^y)$. These works place the problem in a probabilistic framework and tie precise slowly-varying corrections to high moments of $\zeta$.

**Explicit numerical and algorithmic progress (computations and upper bounds for the smallest counterexample):** Because Odlyzko-te Riele's method is non-constructive (it shows a counterexample must exist below some enormous exponential bound), much computational and algorithmic work has gone into two directions: (a) compute $M(n)$ to as high $n$ as practical (giving concrete extrema and improving numeric bounds on $\limsup/\liminf$), and (b) reduce the proven upper bound on the smallest $x$ at which $|M(x)| > \sqrt{x}$. Key contributions:

1. Kotnik & van de Lune (2004) [436] performed numerical experiments on the order of $M(x)$ and gave heuristic evidence and computational data about the local maxima/minima of $M(x)/\sqrt{x}$.

2. Greg Hurst (2018) [437] improved computations of $M(x)$ up to $10^{16}$ and used modern algorithms to push the numerically observed bounds for $\limsup$ and $\liminf$ to roughly $\pm 1.82$ (improving the lower/upper bounds coming from Odlyzko-te Riele methods). Hurst also described improved algorithms for computing $M(x)$ asymptotically faster than naive summation.

3. In comparing to other explicit (rigorous) bounds, the bounds proposed by Cox et. al. (2021) [1214] is weaker than the best known unconditional or conditional bounds in magnitude, but interesting in that it ties the growth of $M(x)$ to $\sqrt{\log(x!)}$ which is much smaller than typical $\sqrt{x}$ scale for large $x$. Its main value is conceptual: exploring novel arithmetic identities and conditional paths toward RH rather than giving the sharpest bounds.

4. More recently (2024-2025) Kim & Nguyen applied advanced lattice-reduction and algorithmic improvements (moving beyond classic LLL to BKZ and modern CVP techniques) to substantially lower the proven upper bound on the smallest counterexample: they report rigorous
reductions of the bound to about $exp(1.96 \times 10^{19})$ (their arXiv / journal work gives the current best rigorous upper bound). These works continue the Odlyzko-te Riele [865] program but exploit decades of progress in lattice algorithms motivated by cryptography.

**Techniques — a brief taxonomy:** The literature splits techniques roughly into these classes:

1. **Explicit formula/complex analysis/Perron's formula:** express $M(x)$ in terms of non-trivial zeta zeros and use zero-location information to bound sums $\sum_{\rho} x^{\rho}/(\rho\zeta'(\rho))$. This is the bridge linking $(M(x))$ to RH and zero statistics.

2. **Zero-density / spacing and moment methods (analytic):** control of zero heights and moments of $\zeta$ feed into conditional upper bounds (Maier, Montgomery, Soundararajan). Soundararajan's argument bounds the frequency of abnormal clustering of zeros and yields the currently best conditional growth bound.

3. **Diophantine/lattice methods (Odlyzko-te Riele and descendants):** translate the existence of large values of the truncated explicit formula into an inhomogeneous simultaneous approximation problem. Then use lattice-basis reduction (LLL, BKZ) and CVP/Aggregation algorithms to produce rigorous upper bounds on the first counterexample. Modern work improves reduction quality and thus the proven exponential bounds.

4. **Heavy computation / algorithmic summation:** faster summation algorithms, GPU methods, blockwise techniques, and practical searches compute $M(x)$ up to very large thresholds and provide empirical data (Kotnik, Hurst, Deléglise-Rivat and others).

Open problems and the current frontier.

1. **Exact order of growth:** Unconditional asymptotic order of $M(x)$ is unknown; conditional heuristics and evidence point to $\sqrt{x}$ times slowly growing iterated-log factors (Gonek/Ng), but no proof exists.

2. **Explicit smallest counterexample:** While Odlyzko-te Riele proved a counterexample exists and later work vastly lowered the proven upper bound for the least counterexample (from astronomically huge values down by many orders of magnitude), an explicit small counterexample is still unknown; current rigorous upper bounds are enormous exponentials (though much reduced by modern lattice work of Kim-Nguyen).

3. **Sharp conditional bounds under RH / GRH:** Soundararajan's conditional bound remains the strongest in the general literature; further progress on zero statistics (moments, spacing) may refine this.

### 19.6.2 Deep Learning Approaches

The explicit formula connecting $M(n)$ to the nontrivial zeros $\rho$ of $\zeta(s)$ is

$$M(n) = \sum_{\rho} \frac{n^{\rho}}{\rho} + O(n). \qquad (19.72)$$

Since the distribution of $\rho$ determines the oscillatory behavior of $M(n)$, deep learning approaches can be employed to analyze these fluctuations. Consider a sequence prediction model for $M(n)$ based on a recurrent neural network (RNN) where the state equation of the hidden layer is

$$h_t = \sigma(Wh_{t-1} + Ux_t + b), \qquad (19.73)$$
where $h_t$ is the hidden state at time step $t$, $x_t$ is the input (previous values of $M(n)$), and $W, U, b$ are learnable parameters. The loss function to optimize is

$$L = \sum_{t=1}^{N} \left| M(n_t) - \hat{M}(n_t) \right|^2. \qquad (19.74)$$

Since $M(n)$ exhibits chaotic fluctuations, it is beneficial to employ **Fourier analysis** in training. The Fourier transform of $M(n)$ is

$$\hat{M}(f) = \int_{-\infty}^{\infty} M(n)e^{-2\pi ifn}dn. \qquad (19.75)$$

A **wavelet transform** decomposition can be applied using basis functions $\psi_{a,b}(t)$ given by

$$\psi_{a,b}(t) = \frac{1}{\sqrt{|a|}} \psi\left(\frac{t-b}{a}\right), \qquad (19.76)$$

which allows us to extract localized frequency behavior in $M(n)$. Since $M(n)$ is driven by the prime factorization structure of $n$, a **graph neural network (GNN)** can be constructed where nodes correspond to prime factors and edges encode multiplicative relationships. The adjacency matrix $A(n)$ is defined such that

$$M(n) = f(A(n)), \qquad (19.77)$$

where $f$ is a deep function approximator trained to estimate $M(n)$. The connection to the **Riemann Hypothesis** suggests that deep reinforcement learning (DRL) can be used to estimate bounds on $M(n)$. Define a reinforcement learning setup where

$$\text{State} = \{M(n_1), M(n_2), \dots, M(n_k)\}, \quad \text{Action} = \{B(n) \text{ adjustments}\}, \quad \text{Reward} = -(|M(n)| - B(n))^2. \qquad (19.78)$$

Using a wavelet-based neural network, the spectral properties of $M(n)$ can be analyzed via the **scaling function** $\phi(t)$ satisfying

$$\phi(t) = \sum_{k} h_{k} \phi(2t - k). \qquad (19.79)$$

A deep convolutional network can extract self-similar patterns in $M(n)$ by training on feature maps defined as

$$F_l(x) = \sigma \left( \sum_i w_i^{(l)} F_{l-1}(x - i) + b_l \right), \qquad (19.80)$$

where $F_l(x)$ represents the activation at layer $l$. The deep learning integration with classical number theory enables empirical conjecture generation, providing insights into bounds on $M(n)$ and prime number distributions.

Here are some experiments that could be conducted using deep learning to analyze the Mertens function $M(n)$. These experiments leverage different machine learning architectures to identify potential patterns, estimate new bounds, or even gain insights into its connection with the Riemann hypothesis.

### 19.6.3 Experiment 1: Sequence Prediction of the Möbius Function

The objective of this experiment is to train a neural network to predict the values of the Möbius function $\mu(n)$ using past values as inputs. The problem is treated as a time-series forecasting task, where the sequence of past Möbius function values is used to estimate the next value. Given the chaotic and pseudo-random behavior of $\mu(n)$, this experiment aims to explore whether neural
networks can uncover hidden regularities or patterns in its structure. The Möbius function $\mu(n)$ is rigorously defined as follows:

$$ \mu(n) = \begin{cases} 1, & \text{if } n \text{ is square-free with an even number of prime factors,} \\ -1, & \text{if } n \text{ is square-free with an odd number of prime factors,} \\ 0, & \text{if } n \text{ is not square-free.} \end{cases} \qquad (19.81) $$

The summatory Möbius function, also known as the Mertens function, is given by:

$$ M(n) = \sum_{k=1}^{n} \mu(k). \qquad (19.82) $$

To construct the dataset, we create input-output pairs such that:

$$ X_n = (\mu(n-10), \mu(n-9), \dots, \mu(n-1)), \quad Y_n = \mu(n). \qquad (19.83) $$

This setup transforms the problem into a sequence prediction task where the goal is to map past values to the next value in the sequence. The probability distribution of $\mu(n)$ is assumed to follow:

$$ P(\mu(n) = -1) = P(\mu(n) = 0) = P(\mu(n) = 1) = \frac{1}{3}, \qquad (19.84) $$

though deviations occur due to arithmetic fluctuations. Neural architectures such as recurrent neural networks (RNNs), long short-term memory (LSTM) networks, and transformer-based models are employed to capture long-range dependencies and potential hidden structures in $\mu(n)$. The optimization function is chosen as:

$$ \mathcal{L} = \sum_{i=1}^{N} (\hat{\mu}(n_i) - \mu(n_i))^2, \qquad (19.85) $$

where $\hat{\mu}(n_i)$ is the predicted value, minimizing the squared error loss. A Fourier spectral analysis is performed by considering the transformed Möbius sequence:

$$ F(\omega) = \sum_{n=1}^{N} \mu(n) e^{-2\pi i \omega n}, \qquad (19.86) $$

which provides insight into hidden periodicities. If the trained model achieves prediction accuracy significantly exceeding random guessing, this would suggest undiscovered regularities in $\mu(n)$. Conversely, failure would reinforce the conjecture that $\mu(n)$ behaves as a pseudo-random function. This experiment probes the learnability of deep arithmetic properties through machine learning techniques and may inspire further investigations into the relationship between deep learning and number theory.

### 19.6.4 Experiment 2: Estimating Growth Bounds on the Mertens function Using Regression Models

The problem of estimating growth bounds on the Mertens function $M(n)$ using regression models necessitates a deeply rigorous mathematical and scientific framework. Given that $M(n)$ is defined as:

$$ M(n) = \sum_{k=1}^{n} \mu(k), $$
where $\mu(k)$ is the Möbius function, a natural question arises regarding the asymptotic behavior and potential bounding of $M(n)$ with respect to known conjectures and heuristic approximations. A key observation is that the classical Mertens conjecture, which suggested that:

$$|M(n)| \le \sqrt{n},$$

was ultimately disproven. However, refined results indicate that:

$$M(n) = O(n^{1/2+\epsilon}) \quad \text{for any } \epsilon > 0 \text{ under the RH.}$$

The goal of this experiment is to determine whether deep learning-based regression models can capture such behavior or even suggest refined empirical bounds beyond known results. To achieve this, we construct a dataset where the input is $n$ and the output is $M(n)$. Given the irregular oscillatory nature of $M(n)$, standard regression models may struggle to converge to meaningful patterns. We therefore consider transformations such as:

$$Y(n) = \log|M(n)|, \quad Y(n) = \frac{M(n)}{n^{1/2}}, \quad Y(n) = \frac{M(n)}{n^{1/2+\epsilon}},$$

and evaluate whether deep learning models can effectively approximate these relationships. The models employed include:

1. Fully connected deep neural networks (DNNs)

2. Convolutional neural networks (CNNs) applied to structured embeddings of $n$

3. Recurrent neural networks (RNNs) incorporating sequential dependencies in arithmetic functions

The training phase involves utilizing known computational results of $M(n)$ for sufficiently large values of $n$, minimizing the loss function:

$$L(\theta) = \sum_{i} (f_{\theta}(n_{i}) - Y(n_{i}))^{2},$$

where $f_{\theta}$ is the deep learning model parameterized by $\theta$. If $\theta^{*}$ represents the optimal parameters, we analyze whether:

$$\sup_{n} |M(n)| \le f_{\theta^{*}}(n)$$

and compare this bound with existing theoretical upper bounds. Furthermore, spectral analysis methods such as Fourier and wavelet transforms are applied to $M(n)$ to decompose its oscillatory structure into dominant frequency components, revealing potential periodic trends. If the neural network can identify new bounding relationships through empirical training, this would suggest avenues for improving theoretical results. The comparison of deep learning-predicted bounds against classical results such as:

$$M(n) = O(n^{1/2+\epsilon}), \quad M(n) = O(n^{1/2}(\log n)^{2}),$$

serves as a benchmark to assess the efficacy of the deep learning approach. Ultimately, the study investigates whether a neural network can establish improved empirical bounds and whether these findings align with rigorous analytic number theory expectations. If successful, this approach might hint at refined hypotheses regarding the asymptotics of $M(n)$, potentially shedding light on deeper connections with the Riemann Hypothesis and related conjectures.
### 19.6.5 Experiment 3: Neural Network Approximation of the Zeta Function Zeros

We aim to explore the potential of deep learning models in approximating the non-trivial zeros of the Riemann zeta function $\zeta(s)$ and their relationship with the behavior of the Mertens function $M(n)$. Given the deep bonds between the M'obius function $\mu(n)$, the Mertens function, and the zeros of $\zeta(s)$, a neural network-based approach could offer new insights into these interconnections. The primary objective is to train a neural network to approximate the locations of non-trivial zeros of $\zeta(s)$ and analyze whether these zeros provide useful information for limiting or predicting the behavior of $M(n)$. Since the critical line conjecture posits that all nontrivial zeros lie on $Re(s) = \frac{1}{2}$, we seek to investigate whether a trained model can detect implicit patterns in these zeros that relate to the growth of $M(n)$.

The first step in this Approach is Data Representation. We have to construct a data set consisting of known non-trivial zeros $\rho_k = \frac{1}{2} + i\gamma_k$ of $\zeta(s)$, where $\gamma_k$ denotes the imaginary part. We then have to generate the corresponding sequences of the values $M(n)$ and $\mu(n)$ to explore the statistical relationships between these number-theoretic functions and the spectral properties of $\zeta(s)$. We then have to explore alternative data representations, such as Fourier-transformed Möbius sequences or wavelet coefficients of $M(n)$, for improved model interpretability.

The second step in this Approach is Neural Network Architectures. We have to Implement **Fourier Transform Neural Networks (FTNNs)** to leverage spectral properties and identify periodicities in $\mu(n)$ and $M(n)$. We then have to utilize **graph neural networks (GNNs)** to model number-theoretic dependencies, representing prime factorization structures as graph embeddings. We then have to experiment with **deep recurrent architectures** (LSTMs, Transformers) to analyze long-range dependencies in sequences of Möbius function values. The third step in this Approach is Training and Evaluation. We have to train the model to map between the zero distribution of $\zeta(s)$ and features of $M(n)$, such as its upper bounds. We then have to validate whether the learned embeddings capture existing theoretical results, such as the classical bound:

$$M(n) = O(n^{1/2+\epsilon}) \quad \text{for any } \epsilon > 0, \text{ assuming the Riemann Hypothesis.} \quad (19.87)$$

Assess model generalization using synthetic test cases and real-world prime counting functions. Some potential insights are:

* **Numerical Evidence for Spectral Connections:** If the model can learn meaningful patterns in the distribution of zeros, this might provide further empirical support for the spectral interpretation of prime number distributions.

* **Predictive Utility:** A model capable of estimating new zero locations could refine our understanding of the error terms in prime number theorems and potentially guide new conjectures in analytic number theory.

* **Deep Learning as a Theoretical Tool:** While deep learning does not offer rigorous mathematical proofs, its ability to approximate highly nonlinear functions could lead to novel heuristic insights, paving the way for new analytical techniques to study the Mertens function and related zeta function properties.

### 19.6.6 Experiment 4: Graph Neural Networks for Prime Factorization Trees

We investigate the potential of Graph Neural Networks (GNNs) to model the factorization structure of integers and its impact on number-theoretic functions such as the Möbius function $\mu(n)$ and the
Mertens function $M(n)$. The prime factorization of an integer inherently forms a tree-like structure, making it a natural candidate for graph-based learning approaches. The function $\mu(n)$ is defined as follows:

$$ \mu(n) = \begin{cases} 1 & \text{if } n \text{ is square-free with an even number of distinct prime factors,} \\ -1 & \text{if } n \text{ is square-free with an odd number of distinct prime factors,} \\ 0 & \text{otherwise.} \end{cases} \qquad (19.88) $$

The cumulative Mertens function is given by:

$$ M(n) = \sum_{k=1}^{n} \mu(k), \qquad (19.89) $$

and understanding its asymptotic behavior is crucial for problems in analytic number theory. The primary objective is to encode the prime factorization of an integer $n$ as a graph and leverage GNNs to predict values of $\mu(n)$ and analyze statistical fluctuations in $M(n)$. Since the properties of $\mu(n)$ are determined entirely by the prime factorization of $n$, we hypothesize that a well-trained GNN could uncover latent structural patterns in these arithmetic functions.

Each number $n$ is represented as a graph $G_n = (V, E)$, where nodes correspond to prime factors, and edges represent multiplicative relationships. Node features such as prime identity, exponent, and relative position within the factorization tree are encoded using adjacency matrices. The eigenvalues of these matrices, given by:

$$ \lambda_i(G_n) = \text{Eigenvalues of } A(G_n), \qquad (19.90) $$

may provide insights into the underlying arithmetic structure. We employ Message Passing Neural Networks (MPNNs), Graph Convolutional Networks (GCNs), and Graph Attention Networks (GATs) to process these number-theoretic structures. The update rule for each node feature representation $h_v^{(t)}$ in an MPNN follows:

$$ h_v^{(t+1)} = \sigma \left( Wh_v^{(t)} + \sum_{u \in \mathcal{N}(v)} W' h_u^{(t)} \right), \qquad (19.91) $$

where $\mathcal{N}(v)$ denotes the neighborhood of $v$, and $\sigma$ is a non-linearity. The model is trained to classify integers according to their Möbius function values and predict bounds on $M(n)$. A key quantity of interest is the supremum bound:

$$ M(n) = O(n^{1/2+\epsilon}), \quad \text{for any } \epsilon > 0, \qquad (19.92) $$

which is compared against the model's empirical predictions. By analyzing learned representations of $\mu(n)$, we may uncover new heuristics for understanding the randomness of arithmetic functions. The spectral properties of factorization graphs could reveal deep connections to the Riemann Hypothesis, given the existing conjecture:

$$ \sum_{n \le x} \mu(n) = O(x^{1/2+\epsilon}), \qquad (19.93) $$

and how machine learning approaches refine these estimates. Furthermore, extensions of this framework could be applied to other multiplicative functions such as Euler's totient function $\phi(n)$ and divisor functions $d(n)$, contributing to the study of number-theoretic graph structures.
### 19.6.7 Experiment 5: Autoencoders for Dimensionality Reduction of Number-Theoretic Functions

We investigate the application of autoencoders in compressing the Möbius function $\mu(n)$ and Mertens function $M(n)$ values into lower-dimensional representations, aiming to extract hidden structural features within their sequences. Given that $\mu(n)$ exhibits seemingly chaotic behavior, a core question is whether its intrinsic structure can be captured in a compressed latent space. The objective is to encode sequences of Möbius function values $\mu(n)$ using autoencoders and analyze whether their latent representations reveal patterns, periodicities, or structural correlations. If meaningful compression is achieved, this could suggest underlying order within number-theoretic sequences. To achieve this, we define an encoder function $E: \mathbb{R}^N \to \mathbb{R}^d$ and a decoder function $D: \mathbb{R}^d \to \mathbb{R}^N$ such that the reconstruction loss is minimized:

$$ \mathcal{L} = \sum_{i=1}^{N} (\mu(n_i) - D(E(\mu(n_i))))^2. \qquad (19.94) $$

We explore different architectures, such as deep feedforward autoencoders, convolutional autoencoders, and variational autoencoders (VAEs) to optimize encoding efficiency. The latent space analysis involves investigating whether the compressed space captures patterns in Möbius function fluctuations by applying principal component analysis (PCA), t-SNE, or UMAP for visualization. We define the correlation function in latent space as:

$$ \rho(d) = \mathbb{E}[E(\mu(n)) \cdot E(\mu(n+d))], \qquad (19.95) $$

to analyze whether there exist long-range dependencies within the compressed representation. In order to quantify the reconstruction error and randomness evaluation, we compute the variance of reconstruction errors to test if $\mu(n)$ is learnable within a low-dimensional manifold:

$$ \sigma^2 = \frac{1}{N} \sum_{i=1}^{N} (\mu(n_i) - D(E(\mu(n_i))))^2. \qquad (19.96) $$

This variance is compared with the expected behavior under random models of $\mu(n)$ to evaluate whether compression introduces meaningful structure or preserves apparent randomness.

If the autoencoder extracts meaningful lower-dimensional features, it may indicate the presence of latent regularities governing the fluctuations of $\mu(n)$. If no significant compression is possible, this reinforces the conjecture that $\mu(n)$ behaves in a pseudorandom manner, aligning with classical heuristics in analytic number theory. Extending this method to divisor functions, Euler's totient function, or Liouville's function could provide further insights into the compressibility and structural complexity of number-theoretic sequences. Moreover, if eigenvalue analysis of the autoencoder's weight matrices reveals non-trivial spectral properties, this could suggest hidden symmetries in arithmetic function distributions. By leveraging deep learning techniques, this experiment aims to explore whether modern neural architectures can uncover patterns in one of the most fundamental yet enigmatic functions of number theory.
# 20 Advanced Architectures

Advanced architectures in deep learning extend beyond traditional feedforward and recurrent networks by incorporating specialized structures for capturing complex patterns in high-dimensional data. Transformers represent a paradigm shift through their attention mechanisms, which dynamically weight input features without relying on recurrence or convolution. The core operation is scaled dot-product attention, where queries $Q \in \mathbb{R}^{n \times d_k}$, keys $K \in \mathbb{R}^{m \times d_k}$, and values $V \in \mathbb{R}^{m \times d_v}$ interact as:

$$ \text{Attention}(Q, K, V) = \text{softmax}\left(\frac{QK^T}{\sqrt{d_k}}\right)V, $$

with multi-head attention concatenating $h$ parallel attention heads to capture diverse feature relationships. The transformer's encoder-decoder architecture stacks $N$ identical layers, each comprising residual connections around attention and position-wise feedforward networks:

$$ \text{FFN}(\mathbf{x}) = \mathbf{W}_2\text{ReLU}(\mathbf{W}_1\mathbf{x} + \mathbf{b}_1) + \mathbf{b}_2, $$

where $\mathbf{W}_1 \in \mathbb{R}^{d_{ff} \times d}$, $\mathbf{W}_2 \in \mathbb{R}^{d \times d_{ff}}$. Positional encodings inject sequential information through sinusoidal functions:

$$ PE_{(pos,2i)} = \sin(pos/10000^{2i/d}), \quad PE_{(pos,2i+1)} = \cos(pos/10000^{2i/d}), $$

enabling permutation-sensitive processing despite attention's inherent permutation invariance.

Graph Neural Networks (GNNs) operate on non-Euclidean data by propagating node features through graph edges. The message-passing framework generalizes convolutional operations for irregular domains, where node embeddings $\mathbf{h}_v^{(l)}$ at layer $l$ aggregate neighbor information:

$$ \mathbf{h}_v^{(l)} = \phi^{(l)} \left( \mathbf{h}_v^{(l-1)}, \bigoplus_{u \in \mathcal{N}(v)} \psi^{(l)}(\mathbf{h}_v^{(l-1)}, \mathbf{h}_u^{(l-1)}, \mathbf{e}_{uv}) \right), $$

with $\phi^{(l)}, \psi^{(l)}$ as learnable functions, $\bigoplus$ a permutation-invariant aggregator (e.g., sum or max), and $\mathbf{e}_{uv}$ edge features. Graph Attention Networks (GATs) refine this via attention coefficients:

$$ \alpha_{uv} = \text{softmax}_u(\text{LeakyReLU}(\mathbf{a}^T[\mathbf{W}\mathbf{h}_u \parallel \mathbf{W}\mathbf{h}_v])) $$

where $\mathbf{a}$ is a learnable attention vector. The expressive power of GNNs is bounded by the Weisfeiler-Lehman graph isomorphism test, with higher-order GNNs overcoming this limitation through ten-sorized neighborhood aggregation.

Diffusion Models learn data distributions via iterative denoising processes. The forward process gradually corrupts data $\mathbf{x}_0$ with Gaussian noise over $T$ steps:

$$ q(\mathbf{x}_t|\mathbf{x}_{t-1}) = \mathcal{N}(\mathbf{x}_t; \sqrt{1-\beta_t}\mathbf{x}_{t-1}, \beta_t\mathbf{I}), $$

while the reverse process learns to invert this corruption:

$$ p_{\theta}(\mathbf{x}_{t-1}|\mathbf{x}_t) = \mathcal{N}(\mathbf{x}_{t-1}; \mu_{\theta}(\mathbf{x}_t, t), \Sigma_{\theta}(\mathbf{x}_t, t)), $$
where $\mu_{\theta}$ predicts the noise component $\epsilon_{\theta}$:

$$ \mu_{\theta}(\mathbf{x}_t, t) = \frac{1}{\sqrt{\alpha_t}} \left( \mathbf{x}_t - \frac{\beta_t}{\sqrt{1 - \bar{\alpha}_t}} \epsilon_{\theta}(\mathbf{x}_t, t) \right). $$

The training objective minimizes the variational lower bound on the negative log-likelihood, equivalent to denoising score matching:

$$ \mathcal{L} = \mathbb{E}_{t, \mathbf{x}_{0,\epsilon}} \left[ \|\epsilon - \epsilon_{\theta}(\sqrt{\bar{\alpha}_t} \mathbf{x}_0 + \sqrt{1 - \bar{\alpha}_t} \epsilon, t)\|^2 \right], $$

with $\bar{\alpha}_t = \prod_{s=1}^{t} (1 - \beta_s)$. These architectures demonstrate how inductive biases, from attention to graph symmetries and diffusion processes, enable deep learning systems to scale across modalities from text to 3D molecular structures.

# 20.1 Transformers and Attention Mechanisms

## 20.1.1 Literature Review of Transformers and Attention Mechanisms

Vaswani et al. [355] introduced the Transformer architecture, replacing recurrent models with a fully attention-based framework for sequence processing. They formulated the self-attention mechanism, mathematically defining query-key-value (QKV) transformations. They proved scalability advantages over RNNs, showing O(1) parallelization benefits and introduced multi-head attention, enabling contextualized embeddings. Nannepagu et al. [356] explored hybrid AI architectures integrating Transformers with deep reinforcement learning (DQN). They developed a theoretical framework for transformer-augmented reinforcement learning and discussed how self-attention refines feature representations for financial time-series prediction. Rose et al. [357] investigated Vision Transformers (ViTs) for cybersecurity applications, examining attention-based anomaly detection. They theoretically compared self-attention with CNN feature extraction and proposed a new loss function for attention weight refinement in cybersecurity detection models. Buehler [358] explored the theoretical interplay between Graph Neural Networks (GNNs) and Transformer architectures. They developed isomorphic self-attention, which preserves graph topological information and introduced graph-structured positional embeddings within Transformer attention. Tabibpour and Madanizadeh [359] investigated Set Transformers as a theoretical extension of Transformers for high-dimensional dynamic systems and introduced permutation-invariant self-attention mechanisms to replace standard Transformers in decision-making tasks and theoretically formalized attention mechanisms for non-sequential data. Kim et al. (2024) [324] developed a Transformer-based anomaly detection framework for video surveillance. They formalized a new spatio-temporal self-attention mechanism to detect anomalies in videos and extended standard Transformer architectures to handle high-dimensional video data. Ghosh (2024) [91] presented a clear and rigorous derivation of the analytical solution of the Burgers equation using the Cole-Hopf transformation, a classical nonlinear-to-linear mapping technique. By employing the heat kernel method, it provides an explicit solution framework that connects nonlinear viscous flow equations to the linear heat equation. The work serves as a foundational first part of a two-part series, offering deep analytical insights into the mathematical structure of Burgers-type nonlinear PDEs. Li and Dong [360] examined Transformer-based attention mechanisms for wireless communication networks. They introduced hybrid spatial and temporal attention layers for large-scale MIMO channel estimation and provided a rigorous mathematical proof of attention-based signal recovery. Asefa and Assabie [361] investigated language-specific adaptations of Transformer-based translation models. They introduced attention mechanism regularization for low-resource language translation and analyzed the impact of different positional encoding strategies on translation quality. Liao and Chen [362] applied transformer architectures to deepfake detection, analyzing self-attention mechanisms for facial feature analysis. They theoretically compared CNNs and ViTs for forgery detection and
introduced attention-head dropout to improve robustness against adversarial attacks. Jiang et. al. [363] proposed a novel Transformer-based approach for medical imaging reconstruction. They introduced Spatial and Channel-wise Transformer (SCFormer) for enhanced attention-based feature aggregation and theoretically extended contrastive learning to Transformer encoders.

## 20.1.2 Analysis of Transformers and Attention Mechanisms

The Transformer model is an advanced neural network architecture fundamentally defined by the self-attention mechanism, which enables global context-aware computations on sequential data. The model processes an input sequence represented by

$$ \mathbf{X} \in \mathbb{R}^{n \times d_{\text{model}}}, \qquad (20.1) $$

where $n$ denotes the sequence length and $d_{\text{model}}$ the embedding dimensionality. Each token in this sequence is projected into three learned spaces—queries $\mathbf{Q}$, keys $\mathbf{K}$, and values $\mathbf{V}$—using the trainable matrices $\mathbf{W}^Q$, $\mathbf{W}^K$, and $\mathbf{W}^V$, such that

$$ \mathbf{Q} = \mathbf{X}\mathbf{W}^Q, \quad \mathbf{K} = \mathbf{X}\mathbf{W}^K, \quad \mathbf{V} = \mathbf{X}\mathbf{W}^V, \qquad (20.2) $$

where $\mathbf{W}^Q, \mathbf{W}^K, \mathbf{W}^V \in \mathbb{R}^{d_{\text{model}} \times d_k}$, with $d_k$ being the dimensionality of queries and keys. The pairwise similarity between tokens is determined by the dot product $\mathbf{Q}\mathbf{K}^\top$, scaled by the factor $\frac{1}{\sqrt{d_k}}$ to ensure numerical stability, yielding the raw attention scores:

$$ \mathbf{S} = \frac{\mathbf{Q}\mathbf{K}^\top}{\sqrt{d_k}}, \qquad (20.3) $$

where $\mathbf{S} \in \mathbb{R}^{n \times n}$. These scores are normalized using the softmax function, producing the attention weights $\mathbf{A}$, where

$$ A_{ij} = \frac{\exp(S_{ij})}{\sum_{k=1}^{n} \exp(S_{ik})}, \qquad (20.4) $$

ensuring $\sum_{j=1}^n A_{ij} = 1$. The output of the attention mechanism is computed as a weighted sum of the values:

$$ \mathbf{Z} = \mathbf{A}\mathbf{V}, \qquad (20.5) $$

where $\mathbf{Z} \in \mathbb{R}^{n \times d_v}$, with $d_v$ being the dimensionality of the value vectors. This process can be expressed compactly as

$$ \text{Attention}(\mathbf{Q}, \mathbf{K}, \mathbf{V}) = \text{softmax}\left(\frac{\mathbf{Q}\mathbf{K}^\top}{\sqrt{d_k}}\right) \mathbf{V}. \qquad (20.6) $$

Multi-head attention extends this mechanism by splitting $\mathbf{Q}, \mathbf{K}, \mathbf{V}$ into $h$ distinct heads, each operating in its subspace. For the $i$-th head:

$$ \text{head}_i = \text{Attention}(\mathbf{Q}_i, \mathbf{K}_i, \mathbf{V}_i) \qquad (20.7) $$

where $\mathbf{Q}_i = \mathbf{X}\mathbf{W}_i^Q, \mathbf{K}_i = \mathbf{X}\mathbf{W}_i^K, \mathbf{V}_i = \mathbf{X}\mathbf{W}_i^V$. The outputs of all heads are concatenated and linearly transformed:

$$ \text{MultiHead}(\mathbf{Q}, \mathbf{K}, \mathbf{V}) = \text{Concat}(\text{head}_1, \dots, \text{head}_h)\mathbf{W}^O, \qquad (20.8) $$

where $\mathbf{W}^O \in \mathbb{R}^{hd_v \times d_{\text{model}}}$. This architecture enables the model to capture multiple types of relationships simultaneously. Positional encodings are added to the input embeddings $\mathbf{X}$ to preserve sequence order. These encodings $\mathbf{P} \in \mathbb{R}^{n \times d_{\text{model}}}$ are defined as:

$$ P_{(\text{pos}, 2i)} = \sin\left(\frac{\text{pos}}{10000^{2i/d_{\text{model}}}}\right), \quad P_{(\text{pos}, 2i+1)} = \cos\left(\frac{\text{pos}}{10000^{2i/d_{\text{model}}}}\right), \qquad (20.9) $$
ensuring unique representations for each position pos and dimension index i. The feedforward network (FFN) applies two dense layers with an intermediate ReLU activation:

$$ FFN(\mathbf{z}) = \max(0, \mathbf{z}\mathbf{W}_1 + \mathbf{b}_1)\mathbf{W}_2 + \mathbf{b}_2, \qquad (20.10) $$

where $\mathbf{W}_1 \in \mathbb{R}^{d_{model} \times d_{ff}}$, $\mathbf{W}_2 \in \mathbb{R}^{d_{ff} \times d_{model}}$, and $d_{ff} > d_{model}$. Residual connections and layer normalization are applied throughout to stabilize training, with the output given by

$$ \mathbf{H}_{\text{out}} = \text{LayerNorm}(\mathbf{H}_{\text{in}} + \text{FFN}(\mathbf{H}_{\text{in}})). \qquad (20.11) $$

Training optimizes the cross-entropy loss over the output distribution:

$$ \mathcal{L} = - \sum_{t=1}^{T} \log P(y_t | y_{<t}, \mathbf{x}), \qquad (20.12) $$

where $P(y_t | y_{<t}, \mathbf{x})$ is modeled using the softmax over the logits $\mathbf{z}_t \mathbf{W}_{\text{out}} + \mathbf{b}_{\text{out}}$, with parameters $\mathbf{W}_{\text{out}}, \mathbf{b}_{\text{out}}$. The Transformer achieves a complexity of $O(n^2 d_k)$ per attention layer due to the computation of $\mathbf{QK}^\top$, yet its parallelization capabilities render it more efficient than recurrent networks. This mathematical formalism, coupled with innovations like sparse attention and dynamic programming, has solidified the Transformer as the cornerstone of modern sequence modeling tasks. While this quadratic complexity poses challenges for very long sequences, it allows for greater parallelization compared to RNNs, which require $O(n)$ sequential steps. Furthermore, the memory complexity of $O(n^2)$ for storing attention weights can be mitigated using sparse approximations or hierarchical attention structures. The Transformer architecture's flexibility and effectiveness stem from its ability to handle diverse tasks by appropriately modifying its components. For example, in Vision Transformers (ViTs), the input sequence is formed by flattening image patches, and the positional encodings capture spatial relationships. In contrast, in sequence-to-sequence tasks like translation, the cross-attention mechanism enables the decoder to focus on relevant parts of the encoder's output.

In conclusion, the Transformer represents a paradigm shift in neural network design, replacing recurrence with attention and enabling unprecedented scalability and performance. The rigorous mathematical foundation of attention mechanisms, combined with the architectural innovations of multi-head attention, positional encoding, and feedforward layers, underpins its success across domains.

## 20.2 Generative Adversarial Networks (GANs)

### 20.2.1 Literature Review of Generative Adversarial Networks (GANs)

Goodfellow et. al. [364] in their landmark paper introduced Generative Adversarial Networks (GANs), where a generator and a discriminator compete in a minimax game. They established the theoretical foundation of adversarial learning and developed the mathematical formulation of GANs using game theory. They also introduced non-cooperative minimax optimization in deep learning. Chappidi and Sundaram [365] extended GANs with graph neural networks (GNNs) for complex real-world perception tasks. They theoretically integrated GANs with reinforcement learning for self-improving models and developed dual Q-learning mechanisms that enhance GAN convergence stability. Joni [366] provided a comprehensive theoretical overview of GAN-based generative models for advanced image synthesis. They formalized GAN loss functions and their optimization challenges and introduced progressive growing GANs as a solution for high-resolution image generation. Li et. al. (2024) [319] extended GANs to materials science, optimizing the crystal structure prediction process. They developed a GAN framework for molecular modeling
and demonstrates GANs in scientific simulations beyond computer vision tasks. Sekhavat (2024) [313] analyzed the philosophy and theoretical basis of GANs in artistic image generation. He discussed GANs from a cognitive science perspective and established a link between GAN training and computational aesthetics. Kalaiarasi and Sudharani (2024) [367] examined GAN-based image steganography, optimizing data hiding techniques using adversarial training. They extended the theoretical properties of adversarial training to security applications and demonstrated how GANs can minimize perceptual distortion in data hiding. Arjmandi-Tash and Mansourian (2024) [368] explored GANs in scientific computing, generating realistic photodetector datasets. They demonstrated GANs as a theoretical tool for synthetic data augmentation and formulated a probabilistic approach to GAN loss functions for sensor modeling. Gao (2024) [369] bridged the gap between GANs and Partial Differential Equations (PDEs) in physics-informed learning. He established a theoretical framework for solving PDEs using GAN-based architectures and developed a new loss function combining adversarial and variational methods. Hisama et. al. [370] applied GANs to computational chemistry, generating new alloy catalyst structures. They introduced Wasserstein GANs (WGANs) for molecular design and used GAN-generated latent spaces to predict catalyst activity. Wang and Zhang (2024) [371] proposed an improved GAN framework for medical image segmentation. They developed a novel attention-enhanced GAN for robust segmentation and provided a mathematical formulation for adversarial segmentation loss functions.

## 20.2.2 Analysis of Generative Adversarial Networks (GANs)

Generative Adversarial Networks (GANs) are an intricate mathematical framework designed to model complex probability distributions by leveraging a competitive dynamic between two neural networks, the generator $G$ and the discriminator $D$. These networks are parametrized by weights $\theta_G \in \Theta_G$ and $\theta_D \in \Theta_D$, and their interaction is mathematically formulated as a two-player zero-sum game. The generator $G: \mathbb{R}^d \to \mathbb{R}^n$ maps latent variables $\mathbf{z} \sim p_{\mathbf{z}}(\mathbf{z})$, where $p_{\mathbf{z}}$ is a prior probability distribution (commonly uniform or Gaussian), to a synthetic data sample $\hat{\mathbf{x}} = G(\mathbf{z})$. The discriminator $D: \mathbb{R}^n \to [0, 1]$ assigns a probability score $D(\mathbf{x})$ indicating whether $\mathbf{x}$ originates from the true data distribution $p_{\text{data}}(\mathbf{x})$ or the generated distribution $p_g(\mathbf{x})$, implicitly defined as the pushforward measure of $p_{\mathbf{z}}$ under $G$, i.e., $p_g = G_{\#}p_{\mathbf{z}}$. The optimization problem governing GANs is expressed as

$$ \min_G \max_D V(D, G) = \mathbb{E}_{\mathbf{x} \sim p_{\text{data}}} [\log D(\mathbf{x})] + \mathbb{E}_{\mathbf{z} \sim p_{\mathbf{z}}} [\log(1 - D(G(\mathbf{z})))], \quad (20.13) $$

where $\mathbb{E}$ denotes the expectation operator. This objective seeks to maximize the discriminator's ability to distinguish between real and generated samples while simultaneously minimizing the generator's ability to produce samples distinguishable from real data. For a fixed generator $G$, the optimal discriminator $D^*$ is obtained by maximizing $V(D, G)$, yielding

$$ D^*(\mathbf{x}) = \frac{p_{\text{data}}(\mathbf{x})}{p_{\text{data}}(\mathbf{x}) + p_g(\mathbf{x})}. \quad (20.14) $$

Substituting $D^*$ back into the value function simplifies it to

$$ V(D^*, G) = \mathbb{E}_{\mathbf{x} \sim p_{\text{data}}} \left[ \log \frac{p_{\text{data}}(\mathbf{x})}{p_{\text{data}}(\mathbf{x}) + p_g(\mathbf{x})} \right] + \mathbb{E}_{\mathbf{x} \sim p_g} \left[ \log \frac{p_g(\mathbf{x})}{p_{\text{data}}(\mathbf{x}) + p_g(\mathbf{x})} \right]. \quad (20.15) $$

This expression is equivalent to minimizing the Jensen-Shannon (JS) divergence between $p_{\text{data}}$ and $p_g$, defined as

$$ \text{JS}(p_{\text{data}}\|p_g) = \frac{1}{2}\text{KL}(p_{\text{data}}\|M) + \frac{1}{2}\text{KL}(p_g\|M), \quad (20.16) $$

where $M = \frac{1}{2}(p_{\text{data}} + p_g)$ and $\text{KL}(p\|q) = \int p(\mathbf{x}) \log \frac{p(\mathbf{x})}{q(\mathbf{x})} d\mathbf{x}$ is the Kullback-Leibler divergence. At the Nash equilibrium, $p_g = p_{\text{data}}$, the JS divergence vanishes, and $D^*(\mathbf{x}) = \frac{1}{2}$ for all $\mathbf{x}$. The gradient
updates during training are derived using stochastic gradient descent. For the discriminator, the gradients are given by

$$ \nabla_{\theta_D} V(D, G) = \mathbb{E}_{\mathbf{x} \sim p_{\text{data}}} [\nabla_{\theta_D} \log D(\mathbf{x})] + \mathbb{E}_{\mathbf{z} \sim p_{\mathbf{z}}} [\nabla_{\theta_D} \log(1 - D(G(\mathbf{z})))]. \quad (20.17) $$

Training Generative Adversarial Networks (GANs) involves iterative updates to the parameters $\theta_D$ of the discriminator and $\theta_G$ of the generator. The discriminator's parameters are updated via gradient ascent to maximize the value function $V(D, G)$, while the generator's parameters are updated via gradient descent to minimize the same value function. Denoting the gradients of $D$ and $G$ with respect to their parameters as $\nabla_{\theta_D}$ and $\nabla_{\theta_G}$, the updates are given by:

$$ \theta_D \leftarrow \theta_D + \eta_D \nabla_{\theta_D} [\mathbb{E}_{x \sim p_{\text{data}}} [\log D(x)] + \mathbb{E}_{z \sim p_z} [\log (1 - D(G(z)))]], \quad (20.18) $$

and

$$ \theta_G \leftarrow \theta_G - \eta_G \nabla_{\theta_G} [\mathbb{E}_{z \sim p_z} [\log (1 - D(G(z)))]]. \quad (20.19) $$

In practice, to address issues of vanishing gradients, an alternative loss function for the generator is often used, defined as:

$$ -\mathbb{E}_{z \sim p_z} [\log D(G(z))] . \quad (20.20) $$

This modification ensures stronger gradient signals when the discriminator is performing well, effectively improving the generator's training dynamics. For the generator, the gradients in the original formulation are expressed as

$$ \nabla_{\theta_G} V(D, G) = -\mathbb{E}_{z \sim p_z} [\nabla_{\theta_G} \log(1 - D(G(\mathbf{z})))], \quad (20.21) $$

but due to vanishing gradients when $D(G(\mathbf{z}))$ is near 0, the non-saturating generator loss is preferred:

$$ \mathcal{L}_G = -\mathbb{E}_{z \sim p_z} [\log D(G(\mathbf{z}))]. \quad (20.22) $$

The convergence of GANs is inherently linked to the properties of $D^*(\mathbf{x})$ and the alignment of $p_g$ with $p_{\text{data}}$. However, mode collapse and training instability are frequently observed due to the non-convex nature of the objective functions. Wasserstein GANs (WGANs) address these issues by replacing the JS divergence with the Wasserstein-1 distance, defined as

$$ W(p_{\text{data}}, p_g) = \inf_{\gamma \in \Pi(p_{\text{data}}, p_g)} \mathbb{E}_{(\mathbf{x}, \mathbf{y}) \sim \gamma} [\|\mathbf{x} - \mathbf{y}\|], \quad (20.23) $$

where $\Pi(p_{\text{data}}, p_g)$ is the set of all couplings of $p_{\text{data}}$ and $p_g$. Using Kantorovich-Rubinstein duality, the Wasserstein distance is reformulated as

$$ W(p_{\text{data}}, p_g) = \sup_{\|f\|_L \le 1} \mathbb{E}_{\mathbf{x} \sim p_{\text{data}}} [f(\mathbf{x})] - \mathbb{E}_{\mathbf{x} \sim p_g} [f(\mathbf{x})], \quad (20.24) $$

where $f$ is a 1-Lipschitz function. To enforce the Lipschitz constraint, gradient penalties are applied, ensuring that $\|\nabla f(\mathbf{x})\| \le 1$.

The mathematical framework of GANs integrates elements from game theory, optimization, and information geometry. Their training involves solving a high-dimensional non-convex game, where theoretical guarantees for convergence are challenging due to saddle points and complex interactions between $G$ and $D$. Nevertheless, GANs represent a mathematically elegant paradigm for generative modeling, with ongoing research extending their theoretical and practical capabilities.
## 20.3 Autoencoders and Variational Autoencoders

### 20.3.1 Literature Review of Autoencoders and Variational Autoencoders

Zhang et. al. (2024) [317] explored a theoretical connection between VAEs and rate-distortion theory in image compression. They established a mathematical framework linking probabilistic autoencoding to lossy image compression and introduced hierarchical variational inference for improving generative modeling capacity. Wang and Huang (2025) [318] developed a formal mathematical proof of convergence in over-parameterized VAEs. They established rigorous mathematical limits for training VAEs and introduced Neural Tangent Kernel (NTK) theory to study how VAEs behave under over-parameterization. Li et. al. (2024) [319] extended VAEs to materials science, optimizing crystal structure prediction. They developed a VAE-based molecular modeling framework and demonstrates the role of generative models beyond image-based applications. Huang (2024) [320] reviewed key techniques in VAEs, GANs, and Diffusion Models for image generation. They analyzed probabilistic modeling in VAEs compared to diffusion-based methods and also established a theoretical hierarchy of generative models. Chenebuah (2024) [321] investigated Autoencoders for energy materials simulation and molecular property prediction. They introduced a novel AE-VAE hybrid model for physical simulations and established a theoretical link between Density Functional Theory (DFT) and Autoencoders. Furth et. al. (2024) [322] explored Graph Neural Networks (GNNs) and VAEs for predicting chemical properties. They established theoretical properties of VAEs for graph-based learning and extended Autoencoders to chemical reaction prediction. Gong et. al. [323] investigated Conditional Variational Autoencoders (CVAEs) for material design. They introduced new loss functions for conditional generative modeling and theoretically proved how VAEs can optimize material selection. Kim et. al. [324] uses Transformer-based Autoencoders (AEs) for video anomaly detection. They established theoretical improvements of AEs for time-series anomaly detection and used spatio-temporal Autoencoder embeddings to capture anomalies in videos. Albert et. al. (2024) [325] compared Kernel Learning Embeddings (KLE) and Variational Autoencoders for dimensionality reduction. They introduced VAE-based models for atmospheric modeling and established a mathematical comparison between VAEs and kernel-based models. Sharma et. al. (2024) [326] explored practical applications of Autoencoders in network intrusion detection. They established Autoencoders as robust feature extractors for anomaly detection and provided a formal study of Autoencoder latent space representations.

### 20.3.2 Analysis of Autoencoders and Variational Autoencoders

An Autoencoder (AE) is an unsupervised learning model that attempts to learn a compact representation of the input data $x \in \mathbb{R}^d$ in a lower-dimensional latent space. This model consists of two primary components: an encoder function $f_{\theta_e}$ and a decoder function $f_{\theta_d}$. The encoder $f_{\theta_e}: \mathbb{R}^d \to \mathbb{R}^l$ maps the input data $x$ to a latent code $z$, where $l \ll d$, representing the compressed information. The decoder $f_{\theta_d}: \mathbb{R}^l \to \mathbb{R}^d$ then reconstructs the input from this latent code, producing an approximation $\hat{x}$. The loss function typically used to train the autoencoder is the reconstruction loss, often formulated as the Mean Squared Error (MSE):

$$ \mathcal{L}(\mathbf{x}, \hat{\mathbf{x}}) = \|\mathbf{x} - \hat{\mathbf{x}}\|_2^2. \qquad (20.25) $$

The optimization procedure seeks to minimize the reconstruction error over the dataset $D$, assuming a distribution $p(\mathbf{x})$ over the input data $\mathbf{x}$. The objective is to learn the optimal parameters $\theta_e$ and $\theta_d$, by solving the following optimization problem:

$$ \min_{\theta_e, \theta_d} \mathbb{E}_{\mathbf{x} \sim p(\mathbf{x})} \left[ \|\mathbf{x} - f_{\theta_d}(f_{\theta_e}(\mathbf{x}))\|_2^2 \right]. \qquad (20.26) $$
This formulation drives the encoder-decoder architecture towards learning a latent representation that preserves key features of the input data, allowing it to be efficiently reconstructed. The solution to this problem is typically pursued via **stochastic gradient descent (SGD)**, where gradients of the loss with respect to the model parameters are computed and backpropagated through the network. In contrast to the deterministic autoencoder, the **Variational Autoencoder (VAE)** introduces a probabilistic model to better capture the distribution of the latent variables. A VAE models the data generation process using a latent variable $\mathbf{z} \in \mathbb{R}^l$, and aims to maximize the likelihood of observing the data $\mathbf{x}$ by integrating over all possible latent variables. Specifically, we have the joint distribution:

$$p(\mathbf{x}, \mathbf{z}) = p(\mathbf{x}|\mathbf{z})p(\mathbf{z}), \qquad (20.27)$$

where $p(\mathbf{x}|\mathbf{z})$ is the likelihood of the data given the latent variables, and $p(\mathbf{z})$ is the prior distribution of the latent variables, typically chosen to be a standard Gaussian $\mathcal{N}(\mathbf{z}; 0, I)$. The prior assumption that $p(\mathbf{z}) = \mathcal{N}(0, I)$ simplifies the modeling, as it imposes no particular structure on the latent space, which allows for flexible modeling of the data distribution. The encoder in a VAE outputs a distribution $q_{\theta_e}(\mathbf{z}|\mathbf{x})$ over the latent variables, typically modeled as a multivariate Gaussian with mean $\mu_{\theta_e}(\mathbf{x})$ and variance $\sigma_{\theta_e}(\mathbf{x})$, i.e. $q_{\theta_e}(\mathbf{z}|\mathbf{x}) = \mathcal{N}(\mathbf{z}; \mu_{\theta_e}(\mathbf{x}), \sigma_{\theta_e}^2(\mathbf{x})I)$. The decoder generates the likelihood of the data $\mathbf{x}$ given the latent variable $\mathbf{z}$, expressed as $p_{\theta_d}(\mathbf{x}|\mathbf{z})$, which typically takes the form of a Gaussian distribution for continuous data. A central challenge in VAE training is the **marginal likelihood** $p(\mathbf{x})$, which represents the probability of the observed data. This marginal likelihood is intractable due to the integral over the latent variables:

$$p(\mathbf{x}) = \int p_{\theta_d}(\mathbf{x}|\mathbf{z})p(\mathbf{z}) d\mathbf{z}. \qquad (20.28)$$

To address this, VAE training employs **variational inference**, which approximates the true posterior $p(\mathbf{z}|\mathbf{x})$ with a variational distribution $q_{\theta_e}(\mathbf{z}|\mathbf{x})$. The goal is to optimize the **Evidence Lower Bound (ELBO)**, which is a lower bound on the log-likelihood $\log p(\mathbf{x})$. The ELBO is derived using **Jensen's inequality**:

$$\log p(\mathbf{x}) \geq \mathbb{E}_{q_{\theta_e}(\mathbf{z}|\mathbf{x})} [\log p_{\theta_d}(\mathbf{x}|\mathbf{z})] - \text{KL}(q_{\theta_e}(\mathbf{z}|\mathbf{x}) || p(\mathbf{z})), \qquad (20.29)$$

where the first term is the expected log-likelihood of the data given the latent variables, and the second term is the **Kullback-Leibler (KL) divergence** between the approximate posterior $q_{\theta_e}(\mathbf{z}|\mathbf{x})$ and the prior $p(\mathbf{z})$. The KL divergence acts as a regularizer, penalizing deviations from the prior distribution. The ELBO can then be written as:

$$\mathcal{L}_{\text{VAE}}(\mathbf{x}) = \mathbb{E}_{q_{\theta_e}(\mathbf{z}|\mathbf{x})} [\log p_{\theta_d}(\mathbf{x}|\mathbf{z})] - \text{KL}(q_{\theta_e}(\mathbf{z}|\mathbf{x}) || p(\mathbf{z})). \qquad (20.30)$$

This formulation balances two competing objectives: maximizing the likelihood of reconstructing $\mathbf{x}$ from $\mathbf{z}$, and minimizing the divergence between the posterior $q_{\theta_e}(\mathbf{z}|\mathbf{x})$ and the prior $p(\mathbf{z})$. In order to perform optimization, we need to compute the gradient of the ELBO with respect to the parameters $\theta_e$ and $\theta_d$. However, since sampling from the distribution $q_{\theta_e}(\mathbf{z}|\mathbf{x})$ is non-differentiable, the **reparameterization trick** is applied. This trick allows us to reparameterize the latent variable $\mathbf{z}$ as:

$$\mathbf{z} = \mu_{\theta_e}(\mathbf{x}) + \sigma_{\theta_e}(\mathbf{x}) \cdot \epsilon, \qquad (20.31)$$

where $\epsilon \sim \mathcal{N}(0, I)$ is a standard Gaussian noise vector. This enables the backpropagation of gradients through the latent space and allows the optimization process to proceed via stochastic gradient descent. In practice, the **Monte Carlo method** is used to estimate the expectation in the ELBO. This involves drawing $K$ samples $\mathbf{z}_k$ from the variational posterior $q_{\theta_e}(\mathbf{z}|\mathbf{x})$ and approximating the expectation as:

$$\hat{\mathcal{L}}_{\text{VAE}}(\mathbf{x}) = \frac{1}{K} \sum_{k=1}^{K} \log p_{\theta_d}(\mathbf{x}|\mathbf{z}_k) - \frac{1}{K} \sum_{k=1}^{K} \log q_{\theta_e}(\mathbf{z}_k|\mathbf{x}). \qquad (20.32)$$
This approximation allows for efficient optimization, even when the latent space is high-dimensional and the exact expectation is computationally prohibitive. Thus, the **training process** of a VAE involves the following steps: first, the encoder produces a distribution $q_{\theta_e}(\mathbf{z}|\mathbf{x})$ for each input $\mathbf{x}$; then, latent variables $\mathbf{z}$ are sampled from this distribution; finally, the decoder reconstructs the data $\hat{\mathbf{x}}$ from the latent variable $\mathbf{z}$. The network is trained to maximize the ELBO, which effectively balances the reconstruction loss and the KL divergence term.

In this rigorous exploration, we have presented the mathematical foundations of both autoencoders and variational autoencoders. The core distinction between the two lies in the introduction of a probabilistic framework in the VAE, which leverages variational inference to optimize a tractable lower bound on the marginal likelihood. Through this process, the VAE learns to generate data by sampling from the latent space and reconstructing the input, while maintaining a well-structured latent distribution through regularization by the KL divergence term. The optimization framework for VAEs is grounded in variational inference and the reparameterization trick, enabling gradient-based optimization techniques to efficiently train deep generative models.

# 20.4 Graph neural networks (GNNs)

## 20.4.1 Literature Review of Graph neural networks (GNNs)

Scarselli et. al. (2009) [464] wrote a foundational paper that introduced the concept of Graph Neural Networks (GNNs). It formalized the idea of processing graph-structured data using neural networks, where nodes iteratively update their representations by aggregating information from their neighbors. The paper laid the theoretical groundwork for GNNs, including convergence guarantees and computational frameworks. Kipf and Welling (2017) [465] introduced Graph Convolutional Networks (GCNs), a simplified and highly effective variant of GNNs. It proposed a localized first-order approximation of spectral graph convolutions, making GNNs scalable and practical for large graphs. GCNs became a cornerstone for many subsequent GNN architectures. Hamilton et. al. (2017) [466] introduced GraphSAGE, a framework for inductive representation learning on large graphs. Unlike transductive methods (e.g., GCN), GraphSAGE generates embeddings for unseen nodes by sampling and aggregating features from a node's local neighborhood. It also introduced mean, LSTM, and pooling aggregators, which are widely used in GNNs. Veličković et. al. (2018) [467] proposed Graph Attention Networks (GATs), which use self-attention mechanisms to compute node representations. GATs assign different weights to neighbors during aggregation, allowing the model to focus on more important nodes. This introduced a new paradigm for handling heterogeneous graph structures. Xu et. al. (2019) [468] analyzed the theoretical expressiveness of GNNs, particularly their ability to distinguish graph structures. It introduced the Graph Isomorphism Network (GIN), which is as powerful as the Weisfeiler-Lehman (WL) graph isomorphism test. The work provided a rigorous framework for understanding the limitations and strengths of GNNs. Gilmer et. al. (2017) [469] formalized the Message Passing Neural Network (MPNN) framework, which generalizes many GNN variants. It introduced a unified view of GNNs as iterative message-passing algorithms, where nodes exchange information with their neighbors. The framework has been widely adopted in molecular and chemical graph analysis. Battaglia et. al. (2018) [470] presented the Graph Networks (GN) framework, which generalizes GNNs to handle relational reasoning over structured data. It introduced a block-based architecture for processing entities, relations, and global attributes, making it applicable to a wide range of tasks, including physics simulations and combinatorial optimization. Bruna et. al. (2014) [471] wrote one of the earliest works that proposed spectral graph convolutions, which use the graph Fourier transform to define convolutional operations on graphs. It laid the foundation for spectral-based GNNs, which later inspired spatial-based methods like GCNs. Ying et. al. (2018) [472] demonstrated the practical application of GNNs in large-scale recommender systems. It introduced PinSage, a GNN-based
model that leverages random walks and efficient sampling techniques to handle web-scale graphs. This work highlighted the scalability and real-world impact of GNNs. Sousa and Ghosh (2025) [532] introduce a new generalization of the Riemann functional equation by deriving an integral representation of the Hurwitz zeta function $\zeta(s, b)$ that eliminates the integral term, establishing a novel relation between $\zeta(s, b)$ and the polylogarithm valid for all complex $s \neq 1$ and positive $b$, leveraging the symmetries between $\zeta(s, b)$ and $\zeta(s, -b)$. Zhou et. al. (2020) [473] wrote a comprehensive review paper that summarized the state-of-the-art in GNNs, covering a wide range of methods, applications, and challenges. It provided a taxonomy of GNN architectures, discussed their theoretical foundations, and highlighted open research directions, making it an essential resource for researchers and practitioners.

## 20.4.2 Analysis of Graph neural networks (GNNs)

Graph Neural Networks (GNNs) are a profound and mathematically intricate class of deep learning models specifically designed to handle and process data that is naturally structured as graphs. Unlike traditional neural networks that operate on Euclidean data structures such as vectors, sequences, or grids, GNNs generalize deep learning to non-Euclidean spaces by directly leveraging the underlying graph topology. The mathematical foundation of GNNs is deeply rooted in algebraic graph theory, spectral graph theory, and the principles of geometric deep learning, all of which contribute to the rigorous understanding of how neural networks can be extended to structured relational data. At the core of any graph-based machine learning model lies the mathematical representation of a graph. Formally, a graph $G$ is defined as an ordered pair $G = (V, E)$, where $V$ represents the set of nodes (or vertices), and $E \subseteq V \times V$ represents the set of edges that define the relationships between nodes. The total number of nodes in the graph is denoted by $|V| = N$, while the number of edges is given by $|E|$. The connectivity of the graph is encoded in the adjacency matrix $A \in \mathbb{R}^{N \times N}$, where $A_{ij}$ is nonzero if and only if there exists an edge between nodes $i$ and $j$. The adjacency matrix can be either binary (indicating the mere presence or absence of an edge) or weighted, in which case $A_{ij}$ encodes the strength or affinity of the connection. In addition to graph connectivity, each node $i$ is often associated with a feature vector $\mathbf{x}_i \in \mathbb{R}^d$, and collecting these feature vectors across all nodes forms the node feature matrix $X \in \mathbb{R}^{N \times d}$, where $d$ is the dimensionality of the feature space.

One of the fundamental challenges in extending neural networks to graph domains is the lack of a consistent node ordering, which makes standard operations such as convolutions, pooling, and fully connected layers non-trivial. Unlike images where a fixed spatial structure allows for well-defined convolutional kernels, graphs exhibit arbitrary structure and permutation invariance, meaning that the labels of nodes can be permuted without altering the intrinsic properties of the graph. This necessitates the development of graph-specific neural network architectures that respect the graph topology while maintaining permutation invariance. To facilitate learning on graphs, GNNs employ a neighborhood aggregation or message-passing scheme, wherein each node iteratively gathers information from its neighbors to update its representation. This process can be formulated mathematically using recursive feature propagation rules. Let $H^{(l)} \in \mathbb{R}^{N \times d_l}$ denote the node feature matrix at layer $l$, where each row $H_i^{(l)}$ represents the embedding of node $i$ at that layer. The most fundamental form of feature propagation follows the update equation:

$$H^{(l+1)} = \sigma \left( \tilde{D}^{-\frac{1}{2}} \tilde{A} \tilde{D}^{-\frac{1}{2}} H^{(l)} W^{(l)} \right), \qquad (20.33)$$

where $W^{(l)} \in \mathbb{R}^{d_l \times d_{l+1}}$ is a learnable weight matrix that transforms the feature representation, $\sigma(\cdot)$ is a nonlinear activation function such as ReLU, and $\tilde{A} = A + I$ is the adjacency matrix augmented with self-loops to ensure that each node includes its own features in the aggregation process. The diagonal matrix $\tilde{D}$ is the degree matrix of $\tilde{A}$, defined as $\tilde{D}_{ii} = \sum_j \tilde{A}_{ij}$, which normalizes the feature propagation to avoid scale distortion. The initial node features are represented as $H^{(0)} = X$, where
X is the matrix of initial node features. The weight matrix at layer l is denoted as $W^{(l)} \in \mathbb{R}^{d_l \times d_{l+1}}$, where $d_l$ and $d_{l+1}$ are the dimensions of the input and output feature spaces at layer l, respectively. The weight matrix is trainable. The adjacency matrix A is augmented with self-loops, denoted as $\tilde{A} = A+I$, where I is the identity matrix. The degree matrix $\tilde{D}$ is the diagonal matrix corresponding to the adjacency matrix with self-loops, defined as:

$$ \tilde{D}_{ii} = \sum_{j} \tilde{A}_{ij}, \qquad (20.34) $$

where $\tilde{A}_{ij}$ is the entry in the augmented adjacency matrix. The function $\sigma(\cdot)$ is a nonlinear activation function (such as ReLU) applied element-wise to the node features. This operation ensures that each node aggregates information from its local neighborhood, facilitating feature propagation across the graph. More generally, GNNs can be defined using a message passing scheme, which consists of two key steps. Each node i receives messages from its neighbors $j \in N(i)$. The aggregated message at node i at layer l is computed as:

$$ m_i^{(l)} = \sum_{j \in N(i)} f_m(H_j^{(l)}, H_i^{(l)}, A_{ij}), \qquad (20.35) $$

where $f_m$ is a learnable function that determines how information is aggregated. The node embedding is updated using the function $f_u$, which takes the current node embedding $H_i^{(l)}$ and the aggregated message $m_i^{(l)}$. The updated embedding for node i at layer $l+1$ is given by:

$$ H_i^{(l+1)} = f_u(H_i^{(l)}, m_i^{(l)}), \qquad (20.36) $$

where $f_u$ is another learnable function. A popular choice for the functions $f_m$ and $f_u$ is:

$$ H_i^{(l+1)} = \sigma \left( W^{(l)} \sum_{j \in N(i)} \tilde{A}_{ij} \tilde{D}_{ii}^{-1} H_j^{(l)} \right), \qquad (20.37) $$

where $W^{(l)}$ is the trainable weight matrix, $\tilde{A}_{ij}$ are the entries of the augmented adjacency matrix, and $\tilde{D}_{ii}^{-1}$ is the inverse of the degree matrix. This formulation is permutation invariant, ensuring that the node embeddings do not depend on the order in which neighbors are processed.

A deeper mathematical understanding of GNNs can be obtained by analyzing their connection to spectral graph theory. The Laplacian matrix, central to spectral graph analysis, is defined as

$$ L = D - A \qquad (20.38) $$

where D is the degree matrix. The normalized Laplacian is given by

$$ L_{\text{sym}} = I - D^{-\frac{1}{2}}AD^{-\frac{1}{2}} \qquad (20.39) $$

which possesses an orthonormal eigenbasis. The eigenvalues of the Laplacian encode fundamental properties of the graph, such as connectivity and diffusion characteristics. Spectral methods define graph convolutions in the Fourier domain using the eigen-decomposition

$$ L_{\text{sym}} = U \Lambda U^{\top} \qquad (20.40) $$

where U is the matrix of eigenvectors and Λ is the diagonal matrix of eigenvalues. The graph Fourier transform of a signal x is then given by

$$ \hat{x} = U^{\top} x \qquad (20.41) $$
and graph convolutions are defined as

$$g_{\theta} * x = U g_{\theta}(\Lambda) U^{\top} x \qquad (20.42)$$

However, this formulation is computationally expensive, requiring the full eigen-decomposition of L, motivating approximations such as Chebyshev polynomials and first-order simplifications like those used in Graph Convolutional Networks (GCNs). Beyond GCNs, several other variants of GNNs have been developed to address limitations and enhance expressivity. Graph Attention Networks (GATs) introduce an attention mechanism to dynamically weight the contributions of neighboring nodes using learnable attention coefficients. The attention mechanism is formulated as:

$$\alpha_{ij} = \frac{\exp(\text{LeakyReLU}(a^{\top}[Wh_i \parallel Wh_j]))}{\sum_{k \in \mathcal{N}(i)} \exp(\text{LeakyReLU}(a^{\top}[Wh_i \parallel Wh_k]))} \qquad (20.43)$$

where $\parallel$ denotes concatenation, $a$ is a learnable parameter vector, and attention scores $\alpha_{ij}$ determine the importance of neighbors in updating node features. Another variant, GraphSAGE, employs different aggregation functions (mean, LSTM-based, or pooling-based) to sample and aggregate information from local neighborhoods, ensuring scalability to large graphs. The theoretical expressivity of GNNs is an active area of research, particularly in the context of the Weisfeiler-Lehman graph isomorphism test. The Graph Isomorphism Network (GIN) is designed to match the expressiveness of the 1-dimensional Weisfeiler-Lehman test, using an aggregation function of the form:

$$H_i^{(l+1)} = \text{MLP} \left( (1 + \epsilon) H_i^{(l)} + \sum_{j \in \mathcal{N}(i)} H_j^{(l)} \right), \qquad (20.44)$$

where MLP(·) is a multi-layer perceptron, and $\epsilon$ is a learnable parameter that controls the contribution of self-information. This formulation has been shown to be more powerful in distinguishing graph structures compared to traditional GCNs. Applications of GNNs span multiple domains, ranging from molecular property prediction in chemistry and biology, where molecules are represented as graphs with atoms as nodes and chemical bonds as edges, to recommendation systems that model users and items as bipartite graphs. Other applications include knowledge graph reasoning, social network analysis, and combinatorial optimization problems.

In summary, Graph Neural Networks represent a mathematically rich extension of deep learning to structured relational data. Their foundation in spectral graph theory, algebraic topology, and geometric deep learning provides a rigorous framework for understanding their function and capabilities. Despite their success, open research challenges remain in improving their expressivity, generalization, and computational efficiency, making them an active and evolving field within modern machine learning.

## 20.5 Physics Informed Neural Networks (PINNs)

### 20.5.1 Literature Review of Physics Informed Neural Networks (PINNs)

Raissi et. al. (2019) [474] wrote a seminal paper that introduced the foundational framework of PINNs. It demonstrates how neural networks can be trained to solve both forward and inverse problems for nonlinear PDEs by incorporating physical laws (e.g., conservation laws, boundary conditions) directly into the loss function. The authors show the effectiveness of PINNs in solving high-dimensional PDEs, such as the Navier-Stokes equations, and highlight their ability to handle noisy and sparse data. Karniadakis et. al. (2021) [475] wrote a review article that provided a comprehensive overview of physics-informed machine learning, with a focus on PINNs. It discusses the theoretical foundations, challenges, and applications of PINNs in solving PDEs,
uncertainty quantification, and data-driven modeling. The paper also highlights the integration of PINNs with other machine learning techniques and their potential for multi-scale and multi-physics problems. Lu et. al. (2021) [476] introduced DeepXDE, a Python library for solving differential equations using deep learning, particularly PINNs. The authors provide a detailed explanation of the library's architecture, its flexibility in handling various types of PDEs, and its ability to solve high-dimensional problems. The paper also includes benchmarks and comparisons with traditional numerical methods. Cox et. al. (2022) [75] introduced a recursive generalization of the classical sum-of-divisors function, extending its analytical structure beyond traditional arithmetic formulations. The proposed generalization yields new classes of superabundant and colossally abundant numbers, offering deeper insights into their growth behavior. Furthermore, the work explores potential implications for Gronwall's theorem and its connections to the Riemann Hypothesis, highlighting rich intersections between number theory and analytic functions. Sirignano and Spiliopoulos (2018) [477] introduced the Deep Galerkin Method (DGM), a precursor to PINNs, which uses deep neural networks to approximate solutions to high-dimensional PDEs. The authors demonstrate the method's effectiveness in solving problems in finance and physics, laying the groundwork for later developments in PINNs. Wang et. al. (2021) [478] addressed a key challenge in training PINNs: the imbalance in gradient flow during optimization, which can lead to poor convergence. The authors propose adaptive weighting schemes and novel architectures to mitigate these issues, significantly improving the robustness and accuracy of PINNs. Mishra and Molinaro (2021) [479] provided a rigorous theoretical analysis of the generalization error of PINNs. The authors derive bounds on the error and discuss the conditions under which PINNs can reliably approximate solutions to PDEs. This paper is crucial for understanding the theoretical limitations and guarantees of PINNs. Zhang et. al. (2020) [480] extended PINNs to solve time-dependent stochastic PDEs by learning in modal space. The authors demonstrate how PINNs can handle uncertainty quantification and stochastic processes, making them applicable to problems in fluid dynamics, materials science, and finance. Jin et. al. (2021) [481] focused on applying PINNs to the incompressible Navier-Stokes equations, a challenging class of PDEs in fluid dynamics. The authors introduce NSFnets, a specialized variant of PINNs, and demonstrate their effectiveness in solving complex flow problems, including turbulent flows. Chen et. al. (2020) [482] showcased the application of PINNs to inverse problems in nano-optics and metamaterials. The authors demonstrate how PINNs can infer material properties and design parameters from limited experimental data, highlighting their potential for real-world engineering applications. Ghosh and Jain (2021) [847] prove the equivalence of Munafo's Conjecture concerning the OEIS sequence A094358 and the sequence formed by all square-free products of prime factors of Fermat numbers (OEIS A023394), thereby establishing a connection between this conjecture and the long-standing assertion that all Fermat numbers are square-free. Although not explicitly about PINNs, the early work of Psichogios and Ungar (1992) [483] laid the groundwork for integrating physical principles with neural networks. It introduces the concept of hybrid modeling, where neural networks are combined with domain knowledge, a precursor to the physics-informed approach used in PINNs.

## 20.5.2 Analysis of Physics Informed Neural Networks (PINNs)

Physics Informed Neural Networks (PINNs) are a class of *neural networks* explicitly designed to incorporate *partial differential equations (PDEs)* and *boundary/initial conditions* into their training process. The goal is to find approximate solutions to the PDEs governing physical systems using *neural networks*, while directly embedding the governing *physical laws* (described by PDEs) into the training mechanism. A typical PDE problem is represented as:

$$ \mathcal{L}u(x) = f(x), \quad \text{for } x \in \Omega \qquad (20.45) $$

where:

* $\mathcal{L}$ is a differential operator, for instance, the Laplace operator $\nabla^2$, or the Navier-Stokes operator for fluid dynamics.
* $u(x)$ is the unknown solution we wish to approximate.

* $f(x)$ is a known source term, which could represent external forces or other sources in the system.

* $\Omega$ is the domain in which the equation is valid, such as a bounded region in $\mathbb{R}^n$ (e.g., $\Omega \subseteq \mathbb{R}^3$).

The solution $u(x)$ is approximated by a neural network $\hat{u}(x, \theta)$, where $\theta$ denotes the parameters (weights and biases) of the neural network. A neural network approximates a function $\hat{u}(x)$ as a composition of nonlinear mappings, typically as:

$$ \hat{u}(x, \theta) = f_{\theta}(x) = \sigma(W_k \sigma(W_{k-1} \cdots \sigma(W_1 x + b_1) + b_2) \cdots + b_k) \quad (20.46) $$

where:

* $\sigma$ is a nonlinear activation function, such as ReLU or sigmoid.

* $W_i$ and $b_i$ are the weight matrices and bias vectors of the i-th layer.

* The function $f_{\theta}(x)$ is a feedforward neural network with multiple layers.

Thus, the neural network learns a representation $\hat{u}(x, \theta)$ that approximates the physical solution to the PDE. The accuracy of this approximation depends on the choice of the network architecture, activation function, and the training process. To enforce that the neural network approximates a solution to the PDE, we introduce a physics-informed loss function. This loss function typically consists of two parts:

1. *Data-driven loss term*: This term enforces the agreement between the model predictions and any available data points (boundary or initial conditions).

2. *Physics-driven loss term*: This term enforces the satisfaction of the governing PDE at collocation points within the domain $\Omega$.

The data-driven component aims to minimize the discrepancy between the predicted solution and the observed values at certain data points. For a set of training data $\{x_i, u_i\}$, the data-driven loss is given by:

$$ \mathcal{L}_{\text{data}} = \sum_{i=1}^{N} |\hat{u}(x_i, \theta) - u_i|^2 \quad (20.47) $$

where $\hat{u}(x_i, \theta)$ is the predicted value of $u(x)$ at $x_i$, and $u_i$ is the observed value.

The physics-driven term ensures that the predicted solution satisfies the PDE. Let $r(x_i)$ represent the PDE residual evaluated at collocation points $x_i \in \Omega$. The residual $r(x_i)$ is defined as the difference between the left-hand side and the right-hand side of the PDE:

$$ r(x_i) = \mathcal{L}\hat{u}(x_i, \theta) - f(x_i) \quad (20.48) $$

Here, $\mathcal{L}\hat{u}(x_i, \theta)$ is the differential operator acting on the neural network approximation $\hat{u}(x_i, \theta)$. By applying automatic differentiation (AD), we can compute the required derivatives of $\hat{u}(x_i, \theta)$ with respect to $x$. For instance, in the case of a second-order differential equation, AD will compute:

$$ \frac{\partial^2 \hat{u}(x)}{\partial x^2} \quad (20.49) $$

The physics-driven loss is then defined as:

$$ \mathcal{L}_{\text{physics}} = \sum_{i=1}^{M} r(x_i)^2 \quad (20.50) $$
where $r(x_i)$ represents the residuals at the collocation points $x_i$ distributed throughout the domain $\Omega$. The number of these points $M$ can vary depending on the problem's complexity and dimensionality. The total loss function is a weighted sum of the data-driven and physics-driven terms:

$$ \mathcal{L}_{\text{total}} = \mathcal{L}_{\text{data}} + \lambda \mathcal{L}_{\text{physics}} \qquad (20.51) $$

where $\lambda$ is a hyperparameter controlling the balance between the two loss terms. Minimizing this loss function during training ensures that the neural network learns to approximate the solution $u(x)$ that satisfies both the data and the governing physical laws. A key feature of PINNs is the use of automatic differentiation (AD), which allows us to compute the derivatives of the neural network approximation $\hat{u}(x, \theta)$ with respect to its inputs (i.e., the spatial coordinates in the PDE). The chain rule of differentiation is applied recursively through the layers of the neural network.

For a neural network with the following architecture:

$$ \hat{u}(x) = f_{\theta}(x) = \sigma(W_k \sigma(\cdots \sigma(W_1 x + b_1) \cdots + b_k)) \qquad (20.52) $$

we can apply backpropagation and automatic differentiation to compute the derivatives $\frac{\partial \hat{u}(x)}{\partial x}$, $\frac{\partial^2 \hat{u}(x)}{\partial x^2}$, and higher derivatives required by the PDE. For example, for the Laplace operator in a 1D setting:

$$ \frac{\partial^2 \hat{u}(x)}{\partial x^2} = \sum_{j=1}^{k} W_j \frac{\partial^2 \sigma(\cdot)}{\partial x^2} \qquad (20.53) $$

This automatic differentiation procedure ensures that the PDE residual $r(x_i) = \mathcal{L}\hat{u}(x_i, \theta) - f(x_i)$ is computed efficiently and accurately. The formulation of PINNs extends naturally to higher-dimensional PDEs. In the case of a system of partial differential equations, the operator $\mathcal{L}$ may involve higher-order derivatives in multiple dimensions. For instance, in fluid dynamics, the governing equations might involve the Navier-Stokes equations, which require computing derivatives in 3D space:

$$ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\nabla p + \nu \nabla^2 \mathbf{u} + \mathbf{f} \qquad (20.54) $$

Here, $\mathbf{u}(x, t)$ is the velocity field, $p$ is the pressure field, and $\nu$ is the kinematic viscosity. The neural network architecture in PINNs can be extended to multi-output networks that solve for vector fields, ensuring that all components of the velocity and pressure fields satisfy the corresponding PDEs. For inverse problems, where we aim to infer unknown parameters of the system (e.g., material properties, boundary conditions), PINNs provide a natural framework. The inverse problem is framed as the optimization of the loss function with respect to both the neural network parameters $\theta$ and the unknown physical parameters $\alpha$:

$$ \mathcal{L}_{\text{total}}(\theta, \alpha) = \mathcal{L}_{\text{data}}(\theta, \alpha) + \lambda \mathcal{L}_{\text{physics}}(\theta) \qquad (20.55) $$

Multi-fidelity PINNs involve using data at multiple levels of fidelity (e.g., coarse vs. fine simulations, experimental data vs. high-fidelity models) to improve training efficiency and accuracy.

Physics-Informed Neural Networks (PINNs) provide an elegant and powerful approach to solving PDEs by embedding physical laws directly into the training process. The use of automatic differentiation allows for efficient computation of residuals, and the combined loss function enforces both data-driven and physics-driven constraints. With applications spanning across many domains in engineering, physics, and biology, PINNs represent a significant advancement in the integration of machine learning with scientific computing.
## 20.6 Implementation of the Deep Galerkin Methods (DGM) using the Physics-Informed Neural Networks (PINNs)

Consider the general form of a partial differential equation (PDE) given by:

$$L(u(x)) = f(x), \quad \text{for } x \in \Omega \qquad (20.56)$$

where $L$ is a differential operator (either linear or nonlinear), $u(x)$ is the unknown solution, and $f(x)$ is a given source or forcing term. The domain $\Omega \subset \mathbb{R}^d$ represents the spatial region where the solution $u(x)$ is sought. In the case of nonlinear PDEs, $L$ may involve both $u(x)$ and its derivatives in a nonlinear fashion. Additionally, boundary conditions are specified as:

$$u(x) = g(x), \quad x \in \partial\Omega \qquad (20.57)$$

where $g(x)$ is a prescribed function at the boundary $\partial\Omega$ of the domain. The weak formulation of the PDE is obtained by multiplying both sides of the differential equation by a test function $v(x)$ and integrating over the domain $\Omega$:

$$\int_{\Omega} v(x)L(u(x)) dx = \int_{\Omega} v(x)f(x) dx \qquad (20.58)$$

This weak formulation is valid in spaces of functions that satisfy appropriate regularity conditions, such as Sobolev spaces. The weak formulation transforms the problem into an integral form, making it easier to handle numerically. The Deep Galerkin Method (DGM) is a deep learning-based method for approximating the solution of PDEs. The fundamental idea is to construct a neural network-based approximation $\hat{u}(x; \theta)$ for the unknown function $u(x)$, such that the residual of the PDE (the error in satisfying the equation) is minimized in a Galerkin sense. This means that the neural network is trained to minimize the weak form of the PDE's residuals over the domain. In the case of DGM using Physics-Informed Neural Networks (PINNs), the solution is embedded in the architecture of a neural network, and the physics of the problem is enforced through the loss function. The PINN aims to minimize the residual of the weak formulation of the PDE, incorporating both the differential equation and boundary conditions. The neural network used to approximate the solution $\hat{u}(x; \theta)$ is typically a feedforward neural network with an input $x \in \mathbb{R}^d$ (where $d$ is the dimension of the domain) and output $\hat{u}(x; \theta)$, which represents the predicted solution at $x$. The parameters $\theta$ represent the weights and biases of the network, and the architecture is chosen to be deep enough to capture the complexity of the solution. The neural network can be expressed as:

$$\hat{u}(x; \theta) = \text{NN}(x; \theta) \qquad (20.59)$$

Here, $\text{NN}(x; \theta)$ denotes the neural network function that maps the input $x$ to an output $\hat{u}(x; \theta)$. The network layers can include nonlinear activation functions such as ReLU or tanh to capture complex behavior. The PINN minimizes a loss function that combines the residual of the PDE and the boundary condition enforcement. Let the loss function be:

$$L(\theta) = L_{\text{PDE}}(\theta) + L_{\text{BC}}(\theta) \qquad (20.60)$$

where $L_{\text{PDE}}(\theta)$ represents the loss due to the PDE residual, and $L_{\text{BC}}(\theta)$ enforces the boundary conditions. The PDE residual $L_{\text{PDE}}(\theta)$ is defined as the error in satisfying the PDE at a set of collocation points $\{x_i\}_{i=1}^{N_{\text{coll}}}$ in the domain $\Omega$, where $N_{\text{coll}}$ is the number of collocation points. The residual at a point $x_i$ is given by the difference between the differential operator applied to the predicted solution $\hat{u}(x_i; \theta)$ and the forcing term $f(x_i)$:

$$L_{\text{PDE}}(\theta) = \frac{1}{N_{\text{coll}}} \sum_{i=1}^{N_{\text{coll}}} (L(\hat{u}(x_i; \theta)) - f(x_i))^2 \qquad (20.61)$$
Here, $L(\hat{u}(x_i; \theta))$ is the result of applying the differential operator to the output of the neural network at the collocation point $x_i$. For nonlinear PDEs, the operator $L$ might involve both $u(x)$ and its derivatives, and the derivatives of $\hat{u}(x; \theta)$ are computed using automatic differentiation. The boundary condition loss $L_{BC}(\theta)$ ensures that the neural network's predictions at boundary points $\{x_{b_i}\}_{i=1}^{N_{BC}}$ satisfy the boundary condition $u(x) = g(x)$. This loss is computed as:

$$L_{BC}(\theta) = \frac{1}{N_{BC}} \sum_{i=1}^{N_{BC}} (\hat{u}(x_{b_i}; \theta) - g(x_{b_i}))^2 \quad (20.62)$$

where $x_{b_i}$ are points on the boundary $\partial\Omega$, and $g(x_{b_i})$ is the prescribed boundary value. For the Training the Neural Network, The objective is to minimize the total loss function:

$$\theta^* = \arg\min_{\theta} (L_{PDE}(\theta) + L_{BC}(\theta)) \quad (20.63)$$

This minimization is achieved using gradient-based optimization algorithms, such as Stochastic Gradient Descent (SGD) or Adam. The gradients of the loss function with respect to the parameters $\theta$ are computed using automatic differentiation, which is a powerful technique in modern deep learning frameworks (e.g., TensorFlow, PyTorch). To achieve a solution in the Galerkin sense, we need to minimize the weak residual of the PDE. The weak residual is derived by integrating the product of the residual and a test function $v(x)$ over the domain:

$$R(x) = L(u(x)) - f(x) \quad (20.64)$$

The weak formulation of the PDE in Galerkin methods ensures that the solution minimizes the projection of the residual onto the space of test functions. In the case of PINNs, the network implicitly constructs this weak form by adjusting the network's parameters to minimize the residual at sampled points. For a general linear PDE, this weak formulation can be expressed as:

$$\int_{\Omega} v(x) (L(\hat{u}(x; \theta)) - f(x)) dx = 0 \quad (20.65)$$

The neural network is designed to minimize the residual $R(x)$ in the weak sense, over the points where the loss is computed.

For nonlinear PDEs, such as the Navier-Stokes equations or nonlinear Schrödinger equations, the neural network's ability to approximate complex functions is key. The operator $L(\hat{u}(x))$ may involve terms like $\hat{u}(x)\nabla\hat{u}(x)$ (nonlinear convection terms), and the neural network can model these nonlinearities by introducing appropriate activation functions in the layers (e.g., ReLU, sigmoid, or tanh). For a nonlinear PDE such as the incompressible Navier-Stokes equations:

$$\frac{\partial \hat{u}}{\partial t} + \hat{u} \cdot \nabla \hat{u} = -\nabla p + \nu \nabla^2 \hat{u} + f \quad (20.66)$$

where $\hat{u}$ is the velocity field, $p$ is the pressure, $\nu$ is the kinematic viscosity, and $f$ is the external forcing, the network learns the solution $\hat{u}(x; \theta)$ and $\hat{p}(x; \theta)$, such that:

$$L(\hat{u}(x; \theta), \hat{p}(x; \theta)) = f(x) \quad (20.67)$$

This requires the network to compute the derivatives of $\hat{u}$ and $\hat{p}$ and use them in the residual computation. Collocation points $x_i$ are typically sampled using Monte Carlo methods or Latin hypercube sampling. This allows for efficient exploration of the domain $\Omega$, especially in high-dimensional spaces. Boundary points $x_{b_i}$ are selected to enforce boundary conditions accurately. The training process uses an iterative optimization procedure (e.g., Adam optimizer) to update the neural network parameters $\theta$. The gradients of the loss function are computed using automatic differentiation.
differentiation in deep learning frameworks, ensuring accurate and efficient computation of the derivatives of $\hat{u}(x)$. Convergence is determined by monitoring the reduction in the total loss $L(\theta)$, which should approach zero as the solution is refined. Residuals are monitored for both the PDE and boundary conditions, ensuring that the solution satisfies the PDE and boundary conditions to a high degree of accuracy.

In the Deep Galerkin Method (DGM) using Physics-Informed Neural Networks (PINNs), we construct a neural network to approximate the solution of a PDE in the weak form. The loss function enforces both the PDE residual and boundary conditions, and the network is trained to minimize this loss using gradient-based optimization. The method is highly flexible and can handle both linear and nonlinear PDEs, leveraging the power of neural networks to solve complex differential equations in a scientifically and mathematically rigorous manner. This rigorous framework can be applied to a wide variety of differential equations, from simple linear cases to complex nonlinear systems, and serves as a powerful tool for solving high-dimensional and difficult-to-solve PDEs.
21 Deep Kolmogorov Methods

Deep Kolmogorov methods refer to a class of techniques that leverage neural networks to approximate solutions to high-dimensional partial differential equations (PDEs), particularly the Kolmogorov PDEs that arise in stochastic processes and mathematical finance. The Kolmogorov backward equation describes the evolution of expectations of a stochastic process $X_t \in \mathbb{R}^d$ governed by the stochastic differential equation (SDE):

$$d\mathbf{X}_t = \mu(\mathbf{X}_t, t)dt + \sigma(\mathbf{X}_t, t)d\mathbf{W}_t,$$

where $\mu$ is the drift, $\sigma$ the diffusion matrix, and $\mathbf{W}_t$ a Wiener process. For a terminal condition $g(\mathbf{X}_T)$, the solution $u(\mathbf{x}, t) = \mathbb{E}[g(\mathbf{X}_T)|\mathbf{X}_t = \mathbf{x}]$ satisfies the PDE:

$$\frac{\partial u}{\partial t} + \mu \cdot \nabla u + \frac{1}{2} \mathrm{Tr}(\sigma \sigma^T \nabla^2 u) = 0, \quad u(\mathbf{x}, T) = g(\mathbf{x}).$$

Traditional numerical methods for solving such PDEs suffer from the curse of dimensionality, but deep learning offers a computationally feasible alternative. The key idea is to approximate $u(\mathbf{x}, t)$ with a neural network $u_\theta(\mathbf{x}, t)$, where $\theta$ denotes the trainable parameters. The network is trained by minimizing a loss function that enforces the PDE and boundary conditions. For the Kolmogorov backward equation, the loss combines the PDE residual and terminal condition:

$$\mathcal{L}(\theta) = \mathbb{E}_{\mathbf{x},t} \left[ \left\| \frac{\partial u_\theta}{\partial t} + \mu \cdot \nabla u_\theta + \frac{1}{2} \mathrm{Tr}(\sigma \sigma^T \nabla^2 u_\theta) \right\|^2 \right] + \mathbb{E}_{\mathbf{x}} \left[ \|u_\theta(\mathbf{x}, T) - g(\mathbf{x})\|^2 \right],$$

where expectations are approximated via Monte Carlo sampling. The gradients $\nabla u_\theta$ and $\nabla^2 u_\theta$ are computed using automatic differentiation, enabling efficient optimization.

A related approach involves the Feynman-Kac formula, which represents the solution $u(\mathbf{x}, t)$ as an expectation over sample paths of the SDE. This motivates algorithms where neural networks learn the mapping from initial conditions to solution profiles by simulating stochastic trajectories. For instance, the deep backward stochastic differential equation (BSDE) method approximates the solution and its gradient simultaneously. The BSDE formulation of the Kolmogorov equation is:

$$dY_t = -f(\mathbf{X}_t, Y_t, Z_t)dt + Z_t^T d\mathbf{W}_t, \quad Y_T = g(\mathbf{X}_T),$$

where $Y_t = u(\mathbf{X}_t, t)$ and $Z_t = \sigma^T \nabla u(\mathbf{X}_t, t)$. A neural network approximates $Z_t$ as $Z_\theta(\mathbf{X}_t, t)$, and the loss function enforces the terminal condition and dynamics:

$$\mathcal{L}(\theta) = \mathbb{E} \left[ \|Y_T - g(\mathbf{X}_T)\|^2 + \int_0^T \|Z_t - Z_\theta(\mathbf{X}_t, t)\|^2 dt \right].$$

The training involves simulating $\mathbf{X}_t$ and jointly optimizing $Y_t$ and $Z_\theta$ to satisfy the BSDE. This method is particularly effective for high-dimensional problems, such as option pricing in finance, where the state space can involve hundreds of dimensions.
Another variant employs physics-informed neural networks (PINNs) to solve Kolmogorov-type PDEs by embedding the differential operator directly into the loss function. The network $u_{\theta}(\mathbf{x}, t)$ is trained to minimize:

$$ \mathcal{L}(\theta) = \mathbb{E}_{\mathbf{x},t} \left[ \left\| \frac{\partial u_{\theta}}{\partial t} + \mathcal{L}_{\text{gen}} u_{\theta} \right\|^2 \right] + \mathbb{E}_{\mathbf{x}} \left[ \| u_{\theta}(\mathbf{x}, 0) - h(\mathbf{x}) \|^2 \right] + \mathbb{E}_{\mathbf{x},t} \left[ \| \mathcal{B} u_{\theta}(\mathbf{x}, t) \|^2 \right], $$

where $\mathcal{L}_{\text{gen}}$ is the generator of the stochastic process, $h(\mathbf{x})$ the initial condition, and $\mathcal{B}$ boundary conditions. The strength of deep Kolmogorov methods lies in their ability to handle high-dimensional problems where traditional grid-based methods fail, leveraging the universal approximation properties of neural networks and modern optimization techniques to provide scalable solutions. These methods have been successfully applied to problems in quantitative finance, stochastic control, and molecular dynamics, demonstrating their versatility and robustness.

## 21.1 Literature Review of Deep Kolmogorov Methods

Han and Jentzen (2017) [497] introduced the Deep BSDE (Backward Stochastic Differential Equation) solver, a foundational framework for solving high-dimensional PDEs using deep learning. It demonstrates how neural networks can approximate solutions to parabolic PDEs by reformulating them as stochastic control problems. The authors rigorously prove the convergence of the method and provide numerical experiments for high-dimensional problems, such as the Hamilton-Jacobi-Bellman equation. Beck et. al. (2021) [498] extended the Deep BSDE method to solve Kolmogorov PDEs, which describe the evolution of probability densities for stochastic processes. The authors provide a theoretical analysis of the approximation capabilities of deep neural networks for these equations and demonstrate the method's effectiveness in high-dimensional settings. While not exclusively focused on Kolmogorov methods, the paper by Raissi et. al. (2019) [474] introduces Physics-Informed Neural Networks (PINNs), which have become a cornerstone in deep learning for PDEs. PINNs incorporate physical laws (e.g., PDEs) directly into the loss function, enabling the solution of forward and inverse problems. The framework is applicable to high-dimensional PDEs and has inspired many subsequent works. Han and Jentzen (2018) [499] provided a comprehensive theoretical and empirical analysis of the Deep BSDE method. It highlights the method's ability to overcome the curse of dimensionality and demonstrates its application to high-dimensional nonlinear PDEs, including those arising in finance and physics. Sirignano and Spiliopoulos (2018) [477] proposed the Deep Galerkin Method (DGM), which uses deep neural networks to approximate solutions to PDEs without requiring a mesh. The method is particularly effective for high-dimensional problems and is shown to outperform traditional numerical methods in certain settings. Yu (2018) [501] introduced the Deep Ritz Method, which uses deep learning to solve variational problems associated with elliptic PDEs. The method is closely related to Kolmogorov methods and provides a powerful alternative for high-dimensional problems. Zhang et. al. (2020) [480] extended PINNs to solve time-dependent stochastic PDEs, including Kolmogorov-type equations. The authors propose a modal decomposition approach to improve the efficiency and accuracy of the method in high dimensions. Jentzen et. al. (2021) [500] provided a rigorous mathematical foundation for deep learning-based methods for nonlinear parabolic PDEs. It includes error estimates and convergence proofs, making it a key reference for understanding the theoretical underpinnings of Deep Kolmogorov Methods. Khoo et. al. (2021) [502] explored the use of neural networks to solve parametric PDEs, which are closely related to Kolmogorov equations. The authors provide a unified framework for handling high-dimensional parameter spaces and demonstrate the method's effectiveness in various applications. While not strictly a deep learning method, the paper by Hutzenthaler et. al. (2020) [503] introduced the Multilevel Picard method, which has inspired many deep learning approaches for high-dimensional PDEs. It provides a theoretical framework for approximating solutions to semilinear parabolic PDEs, including Kolmogorov equations.
The Deep Kolmogorov Method (DKM) is a deep learning-based approach to solving high-dimensional partial differential equations (PDEs), particularly those arising from stochastic processes governed by Itô diffusions. The rigorous foundation of DKM is built upon **stochastic analysis, functional analysis, variational principles, and neural network approximation theory**. To fully understand the method, one must rigorously derive the Kolmogorov backward equation, justify its probabilistic representation via Feynman-Kac theory, and establish the error bounds for deep learning approximations within appropriate function spaces. Let us explore these aspects in their maximal mathematical depth.

## 21.2 The Kolmogorov Backward Equation and Its Functional Formulation

Let $X_t$ be a $d$-dimensional Itô diffusion process defined on a complete filtered probability space $(\Omega, \mathcal{F}, \{\mathcal{F}_t\}_{t \ge 0}, \mathbb{P})$, satisfying the stochastic differential equation (SDE):

$$dX_t = \mu(X_t)dt + \sigma(X_t)dW_t, \quad X_0 = x, \qquad (21.1)$$

where $\mu: \mathbb{R}^d \to \mathbb{R}^d$ is the drift function and $\sigma: \mathbb{R}^d \to \mathbb{R}^{d \times d}$ is the diffusion function. We assume that both $\mu$ and $\sigma$ satisfy global Lipschitz continuity conditions:

$$\|\mu(x) - \mu(y)\| \le C\|x - y\|, \quad \|\sigma(x) - \sigma(y)\| \le C\|x - y\|, \quad \forall x, y \in \mathbb{R}^d. \qquad (21.2)$$

These conditions guarantee the existence of a **unique strong solution** $X_t$ to the SDE, satisfying $\mathbb{E}[\sup_{0 \le t \le T} \|X_t\|^2] < \infty$. The **Kolmogorov backward equation** describes the evolution of a function $u(t, x)$, which is defined as the expected value of a terminal function $g(X_T)$ and an integral source term $f(t, X_t)$:

$$u(t, x) = \mathbb{E} \left[ g(X_T) + \int_t^T f(s, X_s) ds \mid X_t = x \right]. \qquad (21.3)$$

This function satisfies the **parabolic PDE**:

$$\frac{\partial u}{\partial t} + \mathcal{L}u = f, \quad u(T, x) = g(x), \qquad (21.4)$$

where the **second-order differential operator** $\mathcal{L}$, known as the **infinitesimal generator** of $X_t$, is given by:

$$\mathcal{L}u = \sum_{i=1}^{d} \mu_i(x) \frac{\partial u}{\partial x_i} + \frac{1}{2} \sum_{i,j=1}^{d} (\sigma \sigma^T)_{ij} \frac{\partial^2 u}{\partial x_i \partial x_j}. \qquad (21.5)$$

This equation is well-posed in function spaces such as **Sobolev spaces $H^k(\mathbb{R}^d)$**, **Hölder spaces $C^{k,\alpha}(\mathbb{R}^d)$**, or **Bochner spaces $L^p(\Omega; H^k(\mathbb{R}^d))$** under standard **parabolic regularity** assumptions.

## 21.3 The Feynman-Kac Representation and Its Justification

To rigorously justify the probabilistic representation of $u(t, x)$, we define the stochastic process:

$$M_t = g(X_T) + \int_t^T f(s, X_s) ds - u(t, X_t). \qquad (21.6)$$
Applying Itô's Lemma to $u(t, X_t)$, we obtain:

$$dM_t = \left( \frac{\partial u}{\partial t} + \mathcal{L}u - f \right) dt + \nabla u^T \sigma dW_t. \quad (21.7)$$

Since $u$ satisfies the PDE, the drift term vanishes, leaving:

$$dM_t = \nabla u^T \sigma dW_t. \quad (21.8)$$

Taking expectations and noting that the stochastic integral has **zero mean**, we conclude that $M_t$ is a **martingale**, which establishes the **Feynman-Kac representation**:

$$u(t, x) = \mathbb{E} \left[ g(X_T) + \int_t^T f(s, X_s) ds \mid X_t = x \right]. \quad (21.9)$$

To prove the above equation, we assume that $X_t$ is a **diffusion process** satisfying the stochastic differential equation (SDE):

$$dX_t = \mu(X_t, t)dt + \sigma(X_t, t)dW_t, \quad X_0 = x. \quad (21.10)$$

Here $W_t$ is a **standard Brownian motion** on a filtered probability space $(\Omega, \mathcal{F}, (\mathcal{F}_t)_{t \ge 0}, \mathbb{P})$. The drift $\mu(x, t)$ and diffusion $\sigma(x, t)$ are assumed to be **Lipschitz continuous** in $x$ and **measurable** in $t$, ensuring **existence and uniqueness** of a strong solution to the SDE. The filtration $(\mathcal{F}_t)$ is the **natural filtration** of $W_t$, satisfying the **usual conditions** (right-continuity and completeness). We consider the backward **parabolic partial differential equation (PDE):**

$$\frac{\partial u}{\partial t} + \mu(x, t) \frac{\partial u}{\partial x} + \frac{1}{2} \sigma^2(x, t) \frac{\partial^2 u}{\partial x^2} = V(x, t)u + f(x, t), \quad (21.11)$$

with **final condition**:

$$u(x, T) = g(x). \quad (21.12)$$

The Feynman-Kac representation states that:

$$u(x, t) = \mathbb{E} \left[ \int_t^T e^{-\int_t^s V(X_r, r) dr} f(X_s, s) ds + e^{-\int_t^T V(X_r, r) dr} g(X_T) \mid X_t = x \right]. \quad (21.13)$$

This provides a **probabilistic representation** of the solution to the PDE. Let's now revisit some prerequisites from Stochastic Calculus and Functional Analysis. For that, we first discuss the existence of the Stochastic Process $X_t$. The existence of $X_t$ follows from the **standard existence and uniqueness theorem** for SDEs when $\mu(x, t)$ and $\sigma(x, t)$ satisfy the **Lipschitz continuity condition**:

$$|\mu(x, t) - \mu(y, t)| + |\sigma(x, t) - \sigma(y, t)| \le L|x - y|. \quad (21.14)$$

Under these conditions, there exists a unique **strong solution** $X_t$ that is **adapted** to $\mathcal{F}_t$. Let's use the Itô's Lemma for Stochastic Processes, for a sufficiently smooth function $\phi(X_t, t)$, Itô's lemma states:

$$d\phi(X_t, t) = \left( \frac{\partial \phi}{\partial t} + \mu \frac{\partial \phi}{\partial x} + \frac{1}{2} \sigma^2 \frac{\partial^2 \phi}{\partial x^2} \right) dt + \sigma \frac{\partial \phi}{\partial x} dW_t. \quad (21.15)$$

This will be **crucial** in proving the Feynman-Kac formula. Now let us prove the Feynman-Kac Formula. The first step is to define the Stochastic Process $Y_s$. Define:

$$Y_s = e^{-\int_t^s V(X_r, r) dr} u(X_s, s). \quad (21.16)$$

Applying Itô's Lemma to $Y_s$, we expand:

$$dY_s = d \left( e^{-\int_t^s V(X_r, r) dr} u(X_s, s) \right). \quad (21.17)$$
Using the product rule for stochastic calculus:

$$dY_s = e^{-\int_t^s V(X_r, r)dr} du(X_s, s) + u(X_s, s)d\left(e^{-\int_t^s V(X_r, r)dr}\right) + d\left[e^{-\int_t^s V(X_r, r)dr}, u(X_s, s)\right]. \quad (21.18)$$

Applying Itô's formula, we get

$$du(X_s, s) = \left( \frac{\partial u}{\partial t} + \mu \frac{\partial u}{\partial x} + \frac{1}{2} \sigma^2 \frac{\partial^2 u}{\partial x^2} \right) ds + \sigma \frac{\partial u}{\partial x} dW_s. \quad (21.19)$$

Differentiating the exponential term:

$$d\left(e^{-\int_t^s V(X_r, r)dr}\right) = -V(X_s, s)e^{-\int_t^s V(X_r, r)dr} ds. \quad (21.20)$$

Thus:

$$dY_s = e^{-\int_t^s V(X_r, r)dr} \left( \frac{\partial u}{\partial t} + \mu \frac{\partial u}{\partial x} + \frac{1}{2} \sigma^2 \frac{\partial^2 u}{\partial x^2} - Vu \right) ds + e^{-\int_t^s V(X_r, r)dr} \sigma \frac{\partial u}{\partial x} dW_s. \quad (21.21)$$

The second step shall be taking the expectation and using the Martingale Property. Define the process:

$$M_s = \int_t^s e^{-\int_t^r V(X_q, q)dq} \sigma \frac{\partial u}{\partial x} dW_r. \quad (21.22)$$

Since $M_s$ is a stochastic integral, it is a martingale with expectation zero:

$$\mathbb{E}[M_T | X_t] = 0. \quad (21.23)$$

Taking expectations on both sides of the equation for $Y_s$:

$$\mathbb{E}[Y_T | X_t] = Y_t + \mathbb{E} \left[ \int_t^T e^{-\int_t^s V(X_r, r)dr} f ds \mid X_t \right]. \quad (21.24)$$

Using the terminal condition $Y_T = e^{-\int_t^T V(X_r, r)dr} g(X_T)$, we obtain:

$$u(x, t) = \mathbb{E} \left[ \int_t^T e^{-\int_t^s V(X_r, r)dr} f(X_s, s)ds + e^{-\int_t^T V(X_r, r)dr} g(X_T) \mid X_t = x \right]. \quad (21.25)$$

## 21.4 Deep Kolmogorov Method: Neural Network Approximation

The Deep Kolmogorov Method approximates the function $u(t, x)$ using a deep neural network $u_{\theta}(t, x)$, parameterized by $\theta$. The loss function is constructed as:

$$\mathcal{L}(\theta) = \mathbb{E} \left[ \sum_{t=0}^{T} \left( u_{\theta}(t, X_t) - g(X_T) - \int_t^T f(s, X_s) ds \right)^2 \right]. \quad (21.26)$$

The parameters $\theta$ are optimized via stochastic gradient descent (SGD):

$$\theta_{n+1} = \theta_n - \eta \nabla_{\theta} \mathcal{L}(\theta_n), \quad (21.27)$$

where $\eta$ is the learning rate. By the universal approximation theorem, a sufficiently deep network with ReLU activation satisfies:

$$\|u - u_{\theta}\|_{L^2} \le C(L^{-1/2}W^{-1/d}), \quad (21.28)$$
where $L$ is the network depth and $W$ is the network width. Let $u : [0, T] \times \mathbb{R}^d \to \mathbb{R}$ be the exact solution to the Kolmogorov backward equation:

$$ \frac{\partial u}{\partial t} + \mathcal{L}u = f, \quad (t, x) \in [0, T] \times \mathbb{R}^d, \qquad (21.29) $$

where $\mathcal{L}$ is the differential operator:

$$ \mathcal{L}u = \sum_{i=1}^{d} b_i(x) \frac{\partial u}{\partial x_i} + \frac{1}{2} \sum_{i,j=1}^{d} a_{ij}(x) \frac{\partial^2 u}{\partial x_i \partial x_j}, \qquad (21.30) $$

with $b_i(x)$ and $a_{ij}(x)$ satisfying smoothness and uniform ellipticity conditions:

$$ \exists \lambda, \Lambda > 0 \quad \text{such that} \quad \lambda |\xi|^2 \le \sum_{i,j=1}^{d} a_{ij}(x) \xi_i \xi_j \le \Lambda |\xi|^2 \quad \forall \xi \in \mathbb{R}^d. \qquad (21.31) $$

We approximate $u(t, x)$ with a neural network function $u_{\theta}(t, x)$ of the form:

$$ u_{\theta}(t, x) = \sum_{j=1}^{M} \beta_j \sigma(W_j x + b_j), \qquad (21.32) $$

where $\sigma : \mathbb{R} \to \mathbb{R}$ is the activation function, $W_j \in \mathbb{R}^d$, $b_j \in \mathbb{R}$, $\beta_j \in \mathbb{R}$ are trainable parameters, $M$ represents the number of neurons. We seek a bound on the **approximation error**:

$$ \|u - u_{\theta}\|_{H^s}. \qquad (21.33) $$

From Sobolev Space Approximation by Deep Neural Networks. We assume $u \in H^s(\mathbb{R}^d)$ with $s > d/2$, so by the **Sobolev embedding theorem**, we obtain:

$$ H^s(\mathbb{R}^d) \hookrightarrow C^{0,\alpha}(\mathbb{R}^d), \quad \alpha = s - d/2. \qquad (21.34) $$

This ensures $u$ is Hölder continuous, which is crucial for pointwise approximation. From the Universal Approximation in Sobolev Norms, **Barron space theorem** and **error estimates in Sobolev norms**, there exists a neural network $u_{\theta}$ such that:

$$ \inf_{\theta} \|u - u_{\theta}\|_{H^s} \le CW^{-s/d}L^{-s/(2d)}. \qquad (21.35) $$

where $W$ is the network width, $L$ is the depth, $C$ depends on the smoothness of $u$. This error bound refines the classical universal approximation theorem by considering **derivatives up to order $s$**.

To find the Neural Network Approximation Error, let us do the Neural Network Approximation of the Kolmogorov Equation. We now examine the **residual error**:

$$ R_{\theta}(t, x) = \frac{\partial u_{\theta}}{\partial t} + \mathcal{L}u_{\theta} - f. \qquad (21.36) $$

From Sobolev estimates, we obtain:

$$ \|R_{\theta}\|_{L^2} \le CW^{-s/d}L^{-s/(2d)}. \qquad (21.37) $$

This follows from the **regularity of solutions to parabolic PDEs**, specifically that:

$$ \|\mathcal{L}u - \mathcal{L}u_{\theta}\|_{L^2} \le C\|u - u_{\theta}\|_{H^s}. \qquad (21.38) $$

Thus, the overall error is:

$$ \|u - u_{\theta}\|_{H^s} \le CW^{-s/d}L^{-s/(2d)}. \qquad (21.39) $$
To find the Asymptotic Rates of Convergence, For large width $W$ and depth $L$, we analyze the asymptotic behavior:

$$ \lim_{W,L \to \infty} \|u - u_{\theta}\|_{H^s} = 0. \qquad (21.40) $$

Moreover, for fixed computational resources, the optimal allocation satisfies:

$$ W \sim L^{d/s}. \qquad (21.41) $$

This achieves the best rate:

$$ \|u - u_{\theta}\|_{H^s} = \mathcal{O}(L^{-s/(2d)}). \qquad (21.42) $$

We have established that the approximation error for deep neural networks in solving the **Kolmogorov backward equation** satisfies the rigorous bound:

$$ \|u - u_{\theta}\|_{H^s} \leq CW^{-s/d}L^{-s/(2d)}, \qquad (21.43) $$

which follows from **Sobolev theory**, **parabolic PDE regularity**, and **universal approximation in higher-order norms**.

Consider the backward Kolmogorov partial differential equation:

$$ \frac{\partial u}{\partial t} + \mathcal{L}u = f, \quad (t, x) \in [0, T] \times \mathbb{R}^d, \qquad (21.44) $$

where the differential operator $\mathcal{L}$ is:

$$ \mathcal{L}u = \sum_{i=1}^{d} b_i(x) \frac{\partial u}{\partial x_i} + \frac{1}{2} \sum_{i,j=1}^{d} a_{ij}(x) \frac{\partial^2 u}{\partial x_i \partial x_j}. \qquad (21.45) $$

By the *Feynman-Kac representation*, the solution is expressed in terms of an expectation over stochastic trajectories:

$$ u(t, x) = \mathbb{E} \left[ g(X_T) + \int_t^T f(s, X_s) ds \mid X_t = x \right]. \qquad (21.46) $$

where $X_s$ follows the *Itô diffusion*:

$$ dX_s = b(X_s)ds + \sigma(X_s)dW_s \qquad (21.47) $$

for a standard Brownian motion $W_s$. We approximate this expectation using *Monte Carlo sampling*. Given $N$ independent samples $X_T^{(i)} \sim p(x, T)$, the *empirical Monte Carlo estimator* is:

$$ u_N(t, x) = \frac{1}{N} \sum_{i=1}^{N} \left[ g(X_T^{(i)}) + \int_t^T f(s, X_s^{(i)}) ds \right]. \qquad (21.48) $$

The *Monte Carlo sampling error* is the deviation:

$$ E_N = u_N(t, x) - u(t, x). \qquad (21.49) $$

For the Measure-Theoretic Representation of Error. Define the probability space:

$$ (\Omega, \mathcal{F}, \mathbb{P}), \qquad (21.50) $$

where $\Omega$ is the sample space of Brownian paths, $\mathcal{F}$ is the filtration generated by $W_s$, $\mathbb{P}$ is the Wiener measure. The random variable $E_N$ is thus defined over this probability space. By the *Law of Large Numbers (LLN)*, we have

$$ \mathbb{P} \left( \lim_{N \to \infty} E_N = 0 \right) = 1. \qquad (21.51) $$
However, for finite $N$, we quantify the error using advanced probability bounds. Regarding the Asymptotic Analysis of Monte Carlo Error, the expectation of the squared error is:

$$ \mathbb{E}[E_N^2] = \frac{1}{N}\mathrm{Var}\left(g(X_T) + \int_t^T f(s, X_s)ds\right). \qquad (21.52) $$

Applying the Central Limit Theorem (CLT), we obtain the asymptotic distribution:

$$ \sqrt{N}E_N \xrightarrow{d} \mathcal{N}(0, \sigma^2), \qquad (21.53) $$

where:

$$ \sigma^2 = \mathrm{Var}\left(g(X_T) + \int_t^T f(s, X_s)ds\right). \qquad (21.54) $$

Thus, the Monte Carlo error satisfies:

$$ E_N = \mathcal{O}_p\left(\frac{1}{\sqrt{N}}\right). \qquad (21.55) $$

We need to find Refined Error Bounds via Concentration Inequalities. To rigorously bound the error, we employ *Hoeffding's inequality*:

$$ \mathbb{P}(|E_N| \ge \epsilon) \le 2 \exp \left(-\frac{2N\epsilon^2}{\sigma^2}\right). \qquad (21.56) $$

For a higher-order bound, we use the *Berry-Esseen theorem*:

$$ \sup_x \left| \mathbb{P} \left( \frac{\sqrt{N}E_N}{\sigma} \le x \right) - \Phi(x) \right| \le \frac{C}{\sqrt{N}}, \qquad (21.57) $$

where $C$ depends on the third moment:

$$ \mathbb{E}\left[\left|g(X_T) + \int_t^T f(s, X_s)ds - u(t, x)\right|^3\right]. \qquad (21.58) $$

From a Functional Analysis Perspective, we need to find Operator Norm Bounds. Define the Monte Carlo estimator as a *linear operator*:

$$ \mathcal{M}_N : L^2(\Omega) \to \mathbb{R}, \qquad (21.59) $$

such that:

$$ \mathcal{M}_N \phi = \frac{1}{N} \sum_{i=1}^{N} \phi(X_T^{(i)}). \qquad (21.60) $$

The error is then the operator norm deviation:

$$ \|\mathcal{M}_N - \mathbb{E}\|_{L^2} = \mathcal{O}\left(\frac{1}{\sqrt{N}}\right). \qquad (21.61) $$

By the *spectral decomposition of the covariance operator*, the error satisfies:

$$ \|E_N\|_{L^2} \le \frac{\lambda_{\max}^{1/2}}{\sqrt{N}}, \qquad (21.62) $$

where $\lambda_{\max}$ is the largest eigenvalue of the covariance matrix. For a more precise error characterization, we use the *Edgeworth Series* for Higher-Order Expansion:

$$ \mathbb{P}\left(\frac{\sqrt{N}E_N}{\sigma} \le x\right) = \Phi(x) + \frac{\rho_3}{6\sqrt{N}}(1-x^2)\phi(x) + \mathcal{O}\left(\frac{1}{N}\right), \qquad (21.63) $$
where $\rho_3$ is the skewness of $g(X_T)+\int_t^T f(s, X_s)ds$, $\phi(x)$ is the standard normal density. We have now mathematically rigorously proved that the *Monte Carlo sampling error* in the Deep Kolmogorov method satisfies:

$$E_N = \mathcal{O}_p\left(\frac{1}{\sqrt{N}}\right), \qquad (21.64)$$

but with *precise higher-order refinements* via *Berry-Esseen theorem* (finite sample error), *Hoeffding's inequality* (concentration bound), *Functional norm bounds* (operator analysis), *Edgeworth expansion* (higher-order moment corrections). Thus, the optimal error decay rate remains $1/\sqrt{N}$, but the prefactors depend on problem-specific variance and moment conditions.

Therefore, the total approximation error consists of two primary components:

1. **Neural Network Approximation Error:**

$$\|u - u_{\theta}\|_{L^2} \le C(L^{-1/2}W^{-1/d}). \qquad (21.65)$$

2. **Monte Carlo Sampling Error:**

$$\mathcal{O}(N^{-1/2}), \qquad (21.66)$$

where $N$ is the number of samples used in SGD.

Combining these estimates, we obtain:

$$\|u - u_{\theta}\|_{L^2} \le C(L^{-1/2}W^{-1/d} + N^{-1/2}). \qquad (21.67)$$

The Deep Kolmogorov Method (DKM) provides a framework for solving high-dimensional PDEs using deep learning, with rigorous theoretical justification from **stochastic calculus**, **functional analysis**, and **neural network theory**.
## 22 Reinforcement Learning

Reinforcement learning (RL) is a computational framework for learning optimal decision-making policies through interaction with an environment, modeled as a Markov decision process (MDP). An MDP is formally defined by the tuple $(S, A, P, R, \gamma)$, where $S$ is the state space, $A$ the action space, $P(s'|s, a)$ the transition probability kernel, $R(s, a, s')$ the reward function, and $\gamma \in [0, 1)$ the discount factor. The agent's behavior is governed by a policy $\pi: S \to P(A)$, mapping states to distributions over actions, with the objective of maximizing the expected cumulative return:

$$J(\pi) = \mathbb{E}_{\tau \sim \pi} \left[ \sum_{t=0}^{\infty} \gamma^t \mathcal{R}(s_t, a_t, s_{t+1}) \right],$$

where $\tau = (s_0, a_0, s_1, \dots)$ denotes trajectories sampled from the environment under policy $\pi$. The value function $V^\pi(s)$ and action-value function $Q^\pi(s, a)$ quantify the expected return from state $s$ (and action $a$):

$$V^\pi(s) = \mathbb{E}_{\tau \sim \pi} \left[ \sum_{k=0}^{\infty} \gamma^k \mathcal{R}(s_{t+k}, a_{t+k}, s_{t+k+1}) | s_t = s \right], \quad Q^\pi(s, a) = \mathbb{E}_{\tau \sim \pi} \left[ \sum_{k=0}^{\infty} \gamma^k \mathcal{R}(s_{t+k}, a_{t+k}, s_{t+k+1}) | s_t = s, a_t \right]$$

These functions satisfy the Bellman equations:

$$V^\pi(s) = \mathbb{E}_{a \sim \pi(\cdot|s), s' \sim \mathcal{P}(\cdot|s, a)} [\mathcal{R}(s, a, s') + \gamma V^\pi(s')], \quad Q^\pi(s, a) = \mathbb{E}_{s' \sim \mathcal{P}(\cdot|s, a)} [\mathcal{R}(s, a, s') + \gamma \mathbb{E}_{a' \sim \pi(\cdot|s')} [Q^\pi(s', a')]]$$

The optimal policy $\pi^*$ achieves $V^*(s) = \sup_\pi V^\pi(s)$ for all $s \in S$, with the corresponding Bellman optimality equations:

$$V^*(s) = \max_{a \in \mathcal{A}} \mathbb{E}_{s' \sim \mathcal{P}(\cdot|s, a)} [\mathcal{R}(s, a, s') + \gamma V^*(s')], \quad Q^*(s, a) = \mathbb{E}_{s' \sim \mathcal{P}(\cdot|s, a)} \left[ \mathcal{R}(s, a, s') + \gamma \max_{a' \in \mathcal{A}} Q^*(s', a') \right].$$

In practice, RL algorithms approximate these quantities due to computational intractability in large state spaces. Value-based methods, such as Q-learning, iteratively update the action-value function:

$$Q(s_t, a_t) \leftarrow Q(s_t, a_t) + \alpha \left( \mathcal{R}(s_t, a_t, s_{t+1}) + \gamma \max_{a'} Q(s_{t+1}, a') - Q(s_t, a_t) \right),$$

where $\alpha$ is the learning rate. Policy gradient methods directly optimize $J(\pi_\theta)$ by ascending the gradient:

$$\nabla_\theta J(\pi_\theta) = \mathbb{E}_{\tau \sim \pi_\theta} \left[ \sum_{t=0}^{\infty} \gamma^t \nabla_\theta \log \pi_\theta(a_t|s_t) Q^{\pi_\theta}(s_t, a_t) \right],$$

using Monte Carlo estimates. Actor-critic algorithms combine these approaches, maintaining both a policy $\pi_\theta$ (actor) and a value function $V_\phi$ (critic), with the critic reducing variance in policy updates by providing baseline estimates:

$$\nabla_\theta J(\pi_\theta) = \mathbb{E}_{\tau \sim \pi_\theta} \left[ \sum_{t=0}^{\infty} \gamma^t \nabla_\theta \log \pi_\theta(a_t|s_t) \left( \mathcal{R}(s_t, a_t, s_{t+1}) + \gamma V_\phi(s_{t+1}) - V_\phi(s_t) \right) \right].$$
A critical challenge in RL is the exploration-exploitation trade-off, addressed by strategies like ε-greedy policies or entropy regularization. In deep RL, function approximation using neural networks enables generalization across states. The Deep Q-Network (DQN) approximates $Q^*$ with a network $Q_\phi$, minimizing the mean squared Bellman error:

$$ \mathcal{L}(\phi) = \mathbb{E}_{(s,a,r,s') \sim \mathcal{D}} \left[ \left( r + \gamma \max_{a'} Q_{\phi}^-(s', a') - Q_{\phi}(s, a) \right)^2 \right], $$

where $\mathcal{D}$ is a replay buffer and $\phi^-$ a target network. Policy optimization in high-dimensional spaces leverages trust region methods, such as Proximal Policy Optimization (PPO), which clips policy updates to ensure stability:

$$ \mathcal{L}^{\text{CLIP}}(\theta) = \mathbb{E}_t \left[ \min \left( \frac{\pi_{\theta}(a_t|s_t)}{\pi_{\theta_{\text{old}}}(a_t|s_t)} \hat{A}_t, \text{clip} \left( \frac{\pi_{\theta}(a_t|s_t)}{\pi_{\theta_{\text{old}}}(a_t|s_t)}, 1-\epsilon, 1+\epsilon \right) \hat{A}_t \right) \right], $$

where $\hat{A}_t$ is the advantage estimate. These methods have enabled RL to scale to complex tasks, from robotic control to game playing, while theoretical convergence guarantees often rely on assumptions about the MDP's ergodicity and the function approximator's capacity.

## 22.1 Literature Review of Reinforcement Learning

Sutton and Barto (2018) [286] [287] (2021) wrote a definitive textbook on reinforcement learning. It covers the fundamental concepts, including Markov decision processes (MDPs), temporal difference learning, policy gradient methods, and function approximation. The second edition expands on deep reinforcement learning, covering advanced algorithms like DDPG, A3C, and PPO. Bertsekas and Tsitsiklis (1996) [288] laid the theoretical foundation for reinforcement learning by introducing neuro-dynamic programming, an extension of dynamic programming methods for decision-making under uncertainty. It rigorously covers approximate dynamic programming, policy iteration, and value function approximation. Kakade (2003) [289] in his thesis formalized the sample complexity of RL, providing theoretical guarantees for how much data is required for an agent to learn optimal policies. It introduces the PAC-RL (Probably Approximately Correct RL) framework, which has significantly influenced how RL algorithms are evaluated. Szepesvári (2010) [290] presented a rigorous yet concise overview of reinforcement learning algorithms, including value iteration, Q-learning, SARSA, function approximation, and policy gradient methods. It provides deep theoretical insights into convergence proofs and performance bounds. Haarnoja et al. (2018) [291] introduced Soft Actor-Critic (SAC), an off-policy deep reinforcement learning algorithm that maximizes expected reward and entropy simultaneously. It provides a strong theoretical framework for handling exploration-exploitation trade-offs in high-dimensional continuous action spaces. Mnih et al. (2015) [292] introduced Deep Q-Networks (DQN), demonstrating how deep learning can be combined with Q-learning to achieve human-level performance in Atari games. The authors address key challenges in reinforcement learning, including function approximation and stability improvements. Konda and Tsitsiklis (2003) [293]. provided a rigorous theoretical analysis of Actor-Critic methods, which combine policy-based and value-based learning. It formally establishes convergence proofs for actor-critic algorithms and introduces the natural gradient method for policy improvement. Levine (2018) [294] introduced a probabilistic inference framework for reinforcement learning, linking RL to Bayesian inference. It provides a theoretical foundation for maximum entropy reinforcement learning, explaining why entropy-regularized objectives lead to better exploration and stability. Mannor et al. (2022) [295] gave one of the most rigorous mathematical treatments of reinforcement learning theory. It covers several topics: PAC guarantees for RL algorithms, Complexity bounds for exploration, Connections between RL and control theory, Convergence rates of popular RL methods. Borkar (2008) [296] rigorously analyzed stochastic approximation methods, which form the theoretical backbone of RL algorithms like TD-learning, Q-learning, and policy gradient
methods. Borkar provides a dynamical systems perspective to convergence analysis, offering deep mathematical insights.

## 22.2 Key Concepts

Reinforcement Learning (RL) is a branch of machine learning that deals with agents making decisions in an environment to maximize cumulative rewards over time. This formalized decision-making process can be described using concepts such as agents, states, actions, and rewards, all of which are mathematically formulated within the framework of a Markov Decision Process (MDP). The following provides an extremely mathematically rigorous discussion of these key concepts. An agent interacts with the environment by taking actions based on the current state of the environment. The goal of the agent is to maximize the expected cumulative reward over time. A policy $\pi$ is a mapping from states to probability distributions over actions. Formally, the policy $\pi$ can be written as:

$$ \pi: \mathcal{S} \to \mathcal{P}(\mathcal{A}), \qquad (22.1) $$

where $\mathcal{S}$ is the state space, $\mathcal{A}$ is the action space, and $\mathcal{P}(\mathcal{A})$ is the set of probability distributions over the actions. The policy can be either deterministic:

$$ \pi(a_t|s_t) = \begin{cases} 1 & \text{if } a_t = \pi(s_t), \\ 0 & \text{otherwise,} \end{cases} \qquad (22.2) $$

where $\pi(s_t)$ is the action chosen in state $s_t$, or stochastic, in which case the policy assigns a probability distribution over actions for each state $s_t$. The goal of reinforcement learning is to find an optimal policy $\pi^*(s_t)$, which maximizes the expected return (cumulative reward) from any initial state. The optimal policy is defined as:

$$ \pi^*(s_t) = \arg\max_{\pi} \mathbb{E} \left[ \sum_{k=0}^{\infty} \gamma^k r_{t+k} | s_t \right], \qquad (22.3) $$

where $\gamma$ is the discount factor that determines the weight of future rewards, and $\mathbb{E}[\cdot]$ represents the expectation under the policy $\pi$. The optimal policy can be derived from the optimal action-value function $Q^*(s_t, a_t)$, which we define in the next section. The state $s_t \in \mathcal{S}$ describes the current situation of the agent at time $t$, encapsulating all relevant information that influences the agent's decision-making process. The state space $\mathcal{S}$ may be either discrete or continuous. The state transitions are governed by a probability distribution $P(s_{t+1}|s_t, a_t)$, which represents the probability of moving from state $s_t$ to state $s_{t+1}$ given action $a_t$. These transitions satisfy the Markov property, meaning the future state depends only on the current state and action, not the history of previous states or actions:

$$ P(s_{t+1}|s_t, a_t) = P(s_{t+1}|s_t, a_t) \quad \forall s_t, s_{t+1} \in \mathcal{S}, a_t \in \mathcal{A}. \qquad (22.4) $$

Additionally, the transition probabilities satisfy the normalization condition:

$$ \sum_{s_{t+1} \in \mathcal{S}} P(s_{t+1}|s_t, a_t) = 1 \quad \forall s_t, a_t. \qquad (22.5) $$

The state distribution $\rho_t(s_t)$ represents the probability of the agent being in state $s_t$ at time $t$. The state distribution evolves over time according to the transition probabilities:

$$ \rho_{t+k}(s_{t+k}) = \sum_{s_t \in \mathcal{S}} P(s_{t+k}|s_t, a_t) \rho_t(s_t), \qquad (22.6) $$

where $\rho_t(s_t)$ is the initial distribution at time $t$, and $\rho_{t+k}(s_{t+k})$ is the distribution at time $t+k$. An action $a_t$ taken at time $t$ by the agent in state $s_t$ leads to a transition to state $s_{t+1}$ and results in a
reward $r_t$. The agent aims to select actions that maximize its long-term reward. The *action-value function* $Q(s_t, a_t)$ quantifies the expected cumulative reward from taking action $a_t$ in state $s_t$ and following the optimal policy thereafter. It is defined as:

$$Q(s_t, a_t) = \mathbb{E} \left[ \sum_{k=0}^{\infty} \gamma^k r_{t+k} \mid s_t, a_t \right]. \qquad (22.7)$$

The optimal action-value function $Q^*(s_t, a_t)$ satisfies the *Bellman Optimality Equation*:

$$Q^*(s_t, a_t) = R(s_t, a_t) + \gamma \sum_{s_{t+1} \in S} P(s_{t+1} | s_t, a_t) \max_{a_{t+1}} Q^*(s_{t+1}, a_{t+1}). \qquad (22.8)$$

This recursive equation provides the foundation for dynamic programming methods such as *value iteration* and *policy iteration*. The optimal policy $\pi^*(s_t)$ is derived by choosing the action that maximizes the action-value function:

$$\pi^*(s_t) = \arg \max_{a_t \in \mathcal{A}} Q^*(s_t, a_t). \qquad (22.9)$$

The *optimal value function* $V^*(s_t)$, representing the expected return from state $s_t$ under the optimal policy, is given by:

$$V^*(s_t) = \max_{a_t \in \mathcal{A}} Q^*(s_t, a_t). \qquad (22.10)$$

The optimal value function satisfies the Bellman equation:

$$V^*(s_t) = \max_{a_t \in \mathcal{A}} \left[ R(s_t, a_t) + \gamma \sum_{s_{t+1} \in S} P(s_{t+1} | s_t, a_t) V^*(s_{t+1}) \right]. \qquad (22.11)$$

The *reward* $r_t$ at time $t$ is a scalar value that represents the immediate benefit (or cost) the agent receives after taking action $a_t$ in state $s_t$. It is a function $R(s_t, a_t)$ mapping state-action pairs to real numbers:

$$r_t = R(s_t, a_t). \qquad (22.12)$$

The agent's objective is to maximize the cumulative reward, which is given by the *total return* from time $t$:

$$G_t = \sum_{k=0}^{\infty} \gamma^k r_{t+k}. \qquad (22.13)$$

The agent seeks to find a policy $\pi$ that maximizes the expected return. The Bellman equation for the expected return is:

$$V^{\pi}(s_t) = R(s_t, \pi(s_t)) + \gamma \sum_{s_{t+1} \in S} P(s_{t+1} | s_t, \pi(s_t)) V^{\pi}(s_{t+1}). \qquad (22.14)$$

This recursive relation helps in solving for the optimal value function. An RL problem is typically modeled as a *Markov Decision Process (MDP)*, which is defined as the tuple:

$$\mathcal{M} = (S, \mathcal{A}, P, R, \gamma), \qquad (22.15)$$

where:

* $S$ is the state space,

* $\mathcal{A}$ is the action space,

* $P(s_{t+1}|s_t, a_t)$ is the state transition probability,
* $R(s_t, a_t)$ is the reward function,

* $\gamma$ is the discount factor.

The agent's goal is to solve the MDP by finding the optimal policy $\pi^*(s_t)$ that maximizes the cumulative expected reward. Reinforcement Learning provides a powerful framework for decision-making in uncertain environments, where the agent seeks to maximize cumulative rewards over time. The core concepts—agents, states, actions, rewards—are formalized mathematically within the structure of a Markov Decision Process, enabling the application of optimization techniques such as dynamic programming, Q-learning, and policy gradient methods to solve complex decision-making problems.

## 22.3 Deep Q-Learning

### 22.3.1 Literature Review of Deep Q-Learning

Alonso and Arias (2025) [372] rigorously explored the mathematical foundations of Q-learning and its convergence properties. The authors analyze viscosity solutions and the Hamilton-Jacobi-Bellman (HJB) equation, demonstrating how Q-learning approximations align with these principles. The work provides new theoretical guarantees for Q-learning under different function approximation settings. Lu et. al. (2024) [373] proposed a factored empirical Bellman operator to mitigate the curse of dimensionality in Deep Q-learning. The authors provide rigorous theoretical analysis on how factorization reduces complexity while preserving optimality. The study improves the scalability of deep reinforcement learning models. Humayoo (2024) [374] extended Temporal Difference (TD) Learning to deep Q-learning using time-scale separation techniques. It introduces a Q(Δ)-learning approach that improves stability and convergence speed in complex environments. Empirical results validate its performance in Atari benchmarks. Jia et. al. (2024) [375] integrated Deep Q-learning (DQL) with Game Theory for anti-jamming strategies in wireless networks. It provides a rigorous theoretical framework on how multi-agent Q-learning can improve resilience against adversarial attacks. The study introduces multi-armed bandit algorithms and their convergence properties. Chai et. al. (2025) [376] provided a mathematical analysis of transfer learning in non-stationary Markov Decision Processes (MDPs). It extends Deep Q-learning to settings where the environment changes over time, establishing error bounds for Q-learning in these domains. Yao and Gong (2024) [377] developed a resilient Deep Q-network (DQN) model for multi-agent systems (MASs) under Byzantine attacks. The work introduces a novel distributed Q-learning approach with provable robustness against adversarial perturbations. Liu et. al. (2025) [378] introduced SGD-TripleQNet, a multi-Q-learning framework that integrates three Deep Q-networks. The authors provide a mathematical foundation and proof of convergence for their model. The paper bridges reinforcement learning with stochastic gradient descent (SGD) optimization. Masood et. al. (2025) [379] merged Deep Q-learning with Game Theory (GT) to optimize energy efficiency in smart agriculture. It proposes a mathematical model for dynamic energy allocation, proving the existence of Nash equilibria in Q-learning-based decision-making environments. Patrick (2024) [380] bridged economic modeling with Deep Q-learning. It formulates dynamic pricing strategies using deep reinforcement learning (DRL) and provides mathematical proofs on how RL adapts to economic shocks. Mimouni and Avrachenkov (2025) [381] introduced a novel Deep Q-learning algorithm that incorporates the Whittle index, a key concept in optimal stopping problems. It proves convergence bounds and applies the model to email recommender systems, demonstrating improved performance over traditional Q-learning methods.
## 22.3.2 Analysis of Deep Q-Learning

Deep Q-Learning (DQL) is an advanced reinforcement learning (RL) technique where the goal is to approximate the optimal action-value function $Q^*(s, a)$ through the use of deep neural networks. In traditional Q-learning, the action-value function $Q(s, a)$ maps a state-action pair to the expected return or cumulative discounted reward from that state-action pair, under the assumption of following an optimal policy. Formally, the Q-function is defined as:

$$Q(s, a) = \mathbb{E} \left[ \sum_{t=0}^{\infty} \gamma^t r_t \mid s_0 = s, a_0 = a \right] \qquad (22.16)$$

where $\gamma \in [0, 1]$ is the discount factor, which determines the weight of future rewards relative to immediate rewards, and $r_t$ is the reward received at time step $t$. The optimal Q-function $Q^*(s, a)$ satisfies the Bellman optimality equation:

$$Q^*(s, a) = \mathbb{E} \left[ r_t + \gamma \max_{a'} Q^*(s_{t+1}, a') \mid s_0 = s, a_0 = a \right] \qquad (22.17)$$

where $s_{t+1}$ is the next state after taking action $a$ in state $s$, and the maximization term represents the optimal future expected reward. This equation represents the recursive structure of the optimal action-value function, where each Q-value is updated based on the reward obtained in the current step and the maximum future reward expected from the next state. The goal is to learn the optimal Q-function through iterative updates, typically using the Temporal Difference (TD) method. In Deep Q-Learning, the Q-function is approximated by a deep neural network, as directly storing Q-values for every state-action pair is computationally infeasible for large state and action spaces.

Let the approximated Q-function be $Q_\theta(s, a)$, where $\theta$ denotes the parameters (weights and biases) of the neural network that approximates the action-value function. The deep Q-network (DQN) aims to learn $Q_\theta(s, a)$ such that it closely approximates $Q^*(s, a)$ over time. The update of the Q-function follows the TD error principle, where the goal is to minimize the difference between the current Q-values and the target Q-values derived from the Bellman equation. The loss function for training the DQN is given by:

$$L(\theta) = \mathbb{E}_{(s_t, a_t, r_t, s_{t+1}) \sim \mathcal{D}} \left[ (y_t - Q_\theta(s_t, a_t))^2 \right] \qquad (22.18)$$

where $\mathcal{D}$ denotes the experience replay buffer containing previous transitions $(s_t, a_t, r_t, s_{t+1})$. The target $y_t$ for the Q-values is defined as:

$$y_t = r_t + \gamma \max_{a'} Q_{\theta^{-}}(s_{t+1}, a') \qquad (22.19)$$

Here, $\theta^{-}$ represents the parameters of the target network, which is a slowly updated copy of the online network parameters $\theta$. The target network $Q_{\theta^{-}}$ is used to generate stable targets for the Q-value updates, and its parameters are updated periodically by copying the parameters from the online network $\theta$ after every $T$ steps. The idea behind this is to stabilize the training by preventing rapid changes in the Q-values due to feedback loops from the Q-network's predictions. The update rule for the network parameters $\theta$ follows the gradient descent method and is expressed as:

$$\nabla_{\theta} L(\theta) = \mathbb{E}_{(s_t, a_t, r_t, s_{t+1}) \sim \mathcal{D}} \left[ (y_t - Q_{\theta}(s_t, a_t)) \nabla_{\theta} Q_{\theta}(s_t, a_t) \right] \qquad (22.20)$$

where $\nabla_{\theta} Q_{\theta}(s_t, a_t)$ is the gradient of the Q-function with respect to the parameters $\theta$, which is computed using backpropagation through the neural network. This gradient is used to update the parameters of the Q-network to minimize the loss function. In reinforcement learning, the agent must balance exploration (trying new actions) and exploitation (selecting actions that maximize the reward). This is often handled by using an epsilon-greedy policy, where the agent selects a
random action with probability $\epsilon$ and the action with the highest Q-value with probability $1 - \epsilon$. The epsilon value is decayed over time to ensure that, as the agent learns, it shifts from exploration to more exploitation. The epsilon-greedy action selection rule is given by:

$$a_t = \begin{cases} \text{random action,} & \text{with probability } \epsilon \\ \arg\max_a Q_{\theta}(s_t, a), & \text{with probability } 1 - \epsilon \end{cases} \qquad (22.21)$$

This policy encourages the agent to explore different actions at the beginning of training and gradually exploit the learned Q-values as training progresses. The decay of $\epsilon$ typically follows an annealing schedule to balance exploration and exploitation effectively. A critical component in stabilizing training in Deep Q-Learning is the use of experience replay. In standard Q-learning, the updates are based on consecutive transitions, which can lead to high correlations between consecutive data points. This correlation can slow down learning or even lead to instability. Experience replay addresses this issue by storing a buffer of past experiences and sampling random mini-batches from this buffer during training. This breaks the correlation between consecutive samples and results in more stable and efficient updates. Mathematically, the loss function for training the network involves random sampling of transitions $(s_t, a_t, r_t, s_{t+1})$ from the experience replay buffer $\mathcal{D}$, and the update to the Q-values is computed using the Bellman error based on the sampled experiences:

$$L(\theta) = \mathbb{E}_{(s_t, a_t, r_t, s_{t+1}) \sim \mathcal{D}} \left[ \left( r_t + \gamma \max_{a'} Q_{\theta-}(s_{t+1}, a') - Q_{\theta}(s_t, a_t) \right)^2 \right] \qquad (22.22)$$

This method ensures that the Q-values are updated in a way that is less sensitive to the order in which experiences are observed, promoting more stable learning dynamics.

Despite its success, the DQL algorithm can still suffer from certain issues such as overestimation bias and instability due to the maximization step in the Bellman equation. Overestimation bias occurs because the maximization operation $\max_{a'} Q_{\theta-}(s_{t+1}, a')$ tends to overestimate the true value, as the Q-values are updated based on the same Q-network. To address this, Double Q-learning was introduced, which uses two separate Q-networks for action selection and value estimation, reducing overestimation bias. In Double Q-learning, the target Q-value is computed using the following equation:

$$y_t = r_t + \gamma Q_{\theta-}(s_{t+1}, \arg\max_{a'} Q_{\theta}(s_{t+1}, a')) \qquad (22.23)$$

This approach helps to mitigate the overestimation problem by decoupling the action selection from the Q-value estimation process. The value of $\arg\max$ is taken from the online network $Q_{\theta}$, while the Q-value for the next state is estimated using the target network $Q_{\theta-}$. Another extension to improve the DQL framework is Dueling Q-Learning, which decomposes the Q-function into two separate components: the state value function $V_{\theta}(s)$ and the advantage function $A_{\theta}(s, a)$. The Q-function is then expressed as:

$$Q_{\theta}(s, a) = V_{\theta}(s) + A_{\theta}(s, a) \qquad (22.24)$$

This decomposition allows the agent to learn the value of a state $V_{\theta}(s)$ independently of the specific actions, thus reducing the number of parameters needed for learning. This is particularly beneficial in environments where many actions have similar expected rewards, as it enables the agent to focus on identifying the value of states rather than overfitting to individual actions.

In conclusion, Deep Q-Learning is an advanced reinforcement learning method that utilizes deep neural networks to approximate the optimal Q-function, enabling agents to handle large state and action spaces. The mathematical formulation of DQL involves minimizing the loss function based on the temporal difference error, utilizing experience replay to stabilize learning, and using target
networks to prevent instability. Extensions such as Double Q-learning and Dueling Q-learning further improve the performance and stability of the algorithm. Despite its remarkable successes, Deep Q-Learning still faces challenges such as overestimation bias and instability, which have been addressed with innovative modifications to the original algorithm.

### 22.3.3 Analysis of Double Q-Learning

Double Q-Learning is a reinforcement learning algorithm designed to mitigate the problem of maximization bias inherent in standard Q-learning. Maximization bias occurs because the same Q-values are used both for selecting and for evaluating an action, leading to an overestimation of action values. This overestimation arises due to the maximization operator in the Q-learning update rule, which systematically prefers noisy overestimates of action values. The core idea of Double Q-Learning is to decouple the selection of the action from its evaluation by maintaining two separate value function estimates, denoted as $Q^A(s, a)$ and $Q^B(s, a)$. These two estimates are updated in a mutually symmetric and complementary fashion.

The fundamental update rule for standard Q-learning is

$$Q(s_t, a_t) \leftarrow Q(s_t, a_t) + \alpha \left[ r_{t+1} + \gamma \max_a Q(s_{t+1}, a) - Q(s_t, a_t) \right] \quad (22.25)$$

where $\alpha$ is the learning rate and $\gamma$ is the discount factor. The term $\max_a Q(s_{t+1}, a)$ is the source of the maximization bias. Double Q-Learning addresses this by randomly choosing, for each update, which of the two Q-functions, $Q^A$ or $Q^B$, to update. With probability 0.5, $Q^A$ is updated using $Q^B$ to determine the value of the next state, and vice-versa. The specific mathematical formulation for the update is as follows. When updating $Q^A$, the action considered optimal in the next state $s_{t+1}$ is selected using the $Q^A$ function itself:

$$a_{t+1}^* = \arg\max_a Q^A(s_{t+1}, a) \quad (22.26)$$

However, the value of this action is evaluated using the other Q-function, $Q^B$. Thus, the update for $Q^A$ becomes

$$Q^A(s_t, a_t) \leftarrow Q^A(s_t, a_t) + \alpha \left[ r_{t+1} + \gamma Q^B(s_{t+1}, a_{t+1}^*) - Q^A(s_t, a_t) \right] \quad (22.27)$$

Symmetrically, when updating $Q^B$, the optimal action for the next state is selected using $Q^B$:

$$a_{t+1}^{**} = \arg\max_a Q^B(s_{t+1}, a) \quad (22.28)$$

and its value is evaluated using $Q^A$, leading to the update

$$Q^B(s_t, a_t) \leftarrow Q^B(s_t, a_t) + \alpha \left[ r_{t+1} + \gamma Q^A(s_{t+1}, a_{t+1}^{**}) - Q^B(s_t, a_t) \right] \quad (22.29)$$

The expected value of the update target in Double Q-Learning is unbiased because $\mathbb{E}[Q^B(s, a^*)]$ is an unbiased estimator of the true value of action $a^*$, provided that the two Q-functions are learned independently and the noise in their estimates is uncorrelated. This can be expressed as

$$\mathbb{E}[\max_a Q(s, a)] \geq \max_a \mathbb{E}[Q(s, a)] \quad (22.30)$$

with equality only holding if there is no estimation variance. By using two separate functions, the maximization is performed over one set of estimates while the value is taken from the other, ensuring that

$$\mathbb{E}[Q^B(s, \arg\max_a Q^A(s, a))] = \mathbb{E}[Q(s, \arg\max_a Q(s, a))] \quad (22.31)$$
without the upward bias. The overall policy derived from these two value functions is typically based on their sum or average. The action selection policy, often $\epsilon$-greedy, uses the combined Q-function

$$Q^{SUM}(s, a) = Q^{A}(s, a) + Q^{B}(s, a) \qquad (22.32)$$

to select actions, so the greedy action is

$$\pi(s) = \arg \max_a (Q^A(s, a) + Q^B(s, a)) \qquad (22.33)$$

This ensures that the policy is informed by all collected experience. The convergence of Double Q-Learning to the optimal policy $\pi^*$ and the optimal value function $Q^*$ can be proven under standard stochastic approximation conditions, demonstrating that

$$\lim_{t \to \infty} Q^A(s, a) = Q^*(s, a) \qquad (22.34)$$

and

$$\lim_{t \to \infty} Q^B(s, a) = Q^*(s, a) \qquad (22.35)$$

for all state-action pairs $(s, a)$ with probability one, assuming all pairs are visited infinitely often. The decoupling of the selection and evaluation effectively transforms the update into a stochastic approximation algorithm for solving the system of Bellman optimality equations,

$$Q^*(s, a) = \sum_{s'} P(s'|s, a) \left[ R(s, a, s') + \gamma \max_{a'} Q^*(s', a') \right] \qquad (22.36)$$

without the maximization bias that plagues the single-estimator approach.

### 22.3.4 Analysis of Dueling Q-Learning

Dueling Q-Learning, or more precisely the Dueling Network Architecture for Q-Learning, is a deep reinforcement learning innovation that provides a more granular and empirically stable value function approximation by structurally decomposing the state-action value function $Q(s, a)$. The core premise is that for many states, it is unnecessary to estimate the value of each action independently. Instead, the architecture decouples the Q-function into two separate streams: a state value function $V(s)$ and an advantage function $A(s, a)$. The state value function $V(s; \theta, \beta)$ estimates the long-term value of being in state $s$, irrespective of any specific action, formally defined as

$$V^{\pi}(s) = \mathbb{E}_{a_t \sim \pi, s_{t+1} \sim P} \left[ \sum_{t=0}^{\infty} \gamma^t r_t | s_0 = s \right] \qquad (22.37)$$

The advantage function $A(s, a; \theta, \alpha)$ quantifies the relative importance of taking action $a$ in state $s$ compared to the average value of the state, defined as

$$A^{\pi}(s, a) = Q^{\pi}(s, a) - V^{\pi}(s) \qquad (22.38)$$

This relationship is fundamental, as it encapsulates the identity

$$Q^{\pi}(s, a) = V^{\pi}(s) + A^{\pi}(s, a) \qquad (22.39)$$

The neural network architecture implementing this concept consists of a shared feature encoder $f(s; \theta)$ parameterized by $\theta$, which processes the input state $s$. This shared representation is then fed into two separate streams of fully connected layers. One stream, parameterized by weights $\beta$, outputs a scalar $V(s; \theta, \beta)$. The other stream, parameterized by weights $\alpha$, outputs a vector
$A(s, a; \theta, \alpha)$ with one entry for each action $a \in \mathcal{A}$. The naive approach to combining these streams would be to simply compute

$$Q(s, a) = V(s; \theta, \beta) + A(s, a; \theta, \alpha) \qquad (22.40)$$

However, this formulation is unidentifiable. Given the Q-values, one cannot uniquely recover $V(s)$ and $A(s, a)$ because for any constant $c$,

$$Q(s, a) = (V(s) - c) + (A(s, a) + c) \qquad (22.41)$$

This lack of identifiability can lead to poor optimization dynamics and performance in practice.

To enforce a stable and identifiable learning process, the Dueling architecture forces the advantage function estimator to have zero advantage for the chosen action. This is achieved by subtracting the mean advantage of all actions in that state from the raw advantages. The module's aggregate Q-value output is thus computed as

$$Q(s, a; \theta, \alpha, \beta) = V(s; \theta, \beta) + \left( A(s, a; \theta, \alpha) - \frac{1}{|\mathcal{A}|} \sum_{a'} A(s, a'; \theta, \alpha) \right) \qquad (22.42)$$

An alternative, more robust formulation uses the maximum advantage instead of the mean:

$$Q(s, a; \theta, \alpha, \beta) = V(s; \theta, \beta) + \left( A(s, a; \theta, \alpha) - \max_{a' \in \mathcal{A}} A(s, a'; \theta, \alpha) \right) \qquad (22.43)$$

This ensures that the estimated advantage of the optimal action is exactly zero, meaning

$$Q(s, a^*; \theta, \alpha, \beta) = V(s; \theta, \beta) \qquad (22.44)$$

for

$$a^* = \arg \max_{a'} Q(s, a'; \theta, \alpha, \beta) \qquad (22.45)$$

which directly mirrors the definition $Q(s, a^*) = V(s)$. The mean operator is generally preferred as it provides smoother gradients and improves stability, as the max operator can be unstable early in training when the estimates are noisy.

The parameters $\theta$, $\alpha$, and $\beta$ are learned end-to-end by minimizing the standard Q-learning loss, such as the Mean Squared Error (MSE) loss against a target network. For a transition $(s_t, a_t, r_t, s_{t+1})$, the loss $\mathcal{L}(\theta, \alpha, \beta)$ is given by

$$\mathcal{L}(\theta, \alpha, \beta) = \mathbb{E}_{(s,a,r,s') \sim \mathcal{D}} \left[ (y^{DQN} - Q(s, a; \theta, \alpha, \beta))^2 \right] \qquad (22.46)$$

where $y^{DQN}$ is the target value. In the context of Deep Q-Networks (DQN), this target is

$$y^{DQN} = r + \gamma \max_{a'} Q(s', a'; \theta^{-}, \alpha^{-}, \beta^{-}) \qquad (22.47)$$

where $\theta^{-}, \alpha^{-}, \beta^{-}$ are the parameters of a target network that are periodically updated. The gradient of this loss with respect to the parameters flows through both the value stream and the advantage stream, allowing the network to learn which states are valuable (high $V(s)$) and which actions are relatively better or worse than the state average (non-zero $A(s, a)$). The key benefit is that the value stream $V(s)$ can be updated for every action taken in a state $s$, leading to more efficient learning and a more robust value estimate, particularly in environments where the choice of action is frequently inconsequential. This architectural inductive bias allows the agent to generalize learning across actions more effectively than a standard Q-network, which must learn the value of each action in each state independently.
### 22.3.5 Analysis of Prioritized Experience Replay

Prioritized Experience Replay is a critical enhancement to the Deep Q-Network (DQN) algorithm that fundamentally alters the sampling distribution from the experience replay buffer $\mathcal{D}$ from a uniform distribution to a non-uniform distribution that favors transitions with a high expected learning progress, as measured by the magnitude of their temporal-difference (TD) error. The standard DQN algorithm minimizes the sequence of loss functions

$$L_i(\theta_i) = \mathbb{E}_{(s,a,r,s') \sim U(\mathcal{D})} \left[ \left( y_i^{\text{DQN}} - Q(s, a; \theta_i) \right)^2 \right] \quad (22.48)$$

where

$$y_i^{\text{DQN}} = r + \gamma \max_{a'} Q(s', a'; \theta_i^{-}) \quad (22.49)$$

and $U(\mathcal{D})$ denotes a uniform distribution over the replay buffer. The core premise of prioritized replay is that not all transitions are equally valuable; some transitions, particularly those with a large TD error

$$\delta = y_i^{\text{DQN}} - Q(s, a; \theta) \quad (22.50)$$

provide a more significant learning signal and should be replayed more frequently. The sampling probability $P(i)$ for a transition $i$ is therefore defined as a function of its priority $p_i$. Two primary schemes for defining $p_i$ are the proportional prioritization, given by

$$p_i = |\delta_i| + \epsilon \quad (22.51)$$

where $\epsilon$ is a small positive constant to ensure non-zero probability, and the rank-based prioritization, given by $p_i = \frac{1}{\text{rank}(i)}$, where $\text{rank}(i)$ is the rank of transition $i$ when the replay buffer is sorted according to $|\delta_i|$. The sampling probability is then formalized as

$$P(i) = \frac{p_i^{\alpha}}{\sum_{k=1}^{|\mathcal{D}|} p_k^{\alpha}} \quad (22.52)$$

where the exponent $\alpha$, with $0 \le \alpha \le 1$, determines the intensity of the prioritization. When $\alpha = 0$, the sampling becomes uniform, and as $\alpha$ approaches 1, it becomes fully greedy with respect to the priorities. This non-uniform sampling introduces a bias into the gradient estimates, as the expectation in the loss function is now drawn from a distribution that depends on the changing parameters $\theta$, violating the i.i.d. assumption central to stochastic gradient descent.

To correct for this sampling bias, the algorithm employs importance-sampling (IS) weights. The ideal IS weight that corrects for the non-uniform sampling is

$$w_i = \left( \frac{1}{|\mathcal{D}|} \cdot \frac{1}{P(i)} \right) \quad (22.53)$$

However, to increase stability and ensure the weights are bounded, this is typically modified by exponentiating by a parameter $\beta$ that is annealed from an initial value $\beta_0$ to 1 over the course of training:

$$w_i = \left( \frac{1}{|\mathcal{D}|} \cdot \frac{1}{P(i)} \right)^{\beta} \quad (22.54)$$

These weights are then normalized by $\frac{1}{\max_j w_j}$ to scale the updates downwards, ensuring the correction is always a contraction. The full prioritized replay loss function thus becomes

$$L_{\text{prioritized}}(\theta) = \mathbb{E}_{i \sim P(i)} \left[ w_i \cdot \left( y_i^{\text{DQN}} - Q(s_i, a_i; \theta) \right)^2 \right] \quad (22.55)$$