is **compact** whenever $q \le \frac{np}{n-p}$. More precisely, this means that if $\{u_k\} \subset W^{1,p}(\Omega)$ is a bounded sequence in the Sobolev norm, i.e., there exists a constant $C > 0$ such that

$$ \|u_k\|_{W^{1,p}(\Omega)} = \|u_k\|_{L^p(\Omega)} + \|\nabla u_k\|_{L^p(\Omega)} \le C, \qquad (1.78) $$

then there exists a **subsequence** $\{u_{k_j}\}$ and a function $u \in L^q(\Omega)$ such that

$$ u_{k_j} \to u \quad \text{strongly in } L^q(\Omega), \qquad (1.79) $$

which means that

$$ \|u_{k_j} - u\|_{L^q(\Omega)} \to 0 \quad \text{as } j \to \infty. \qquad (1.80) $$

To establish this rigorously, we first recall the fact that **bounded sequences in $W^{1,p}(\Omega)$ are weakly precompact**. Since $W^{1,p}(\Omega)$ is a **reflexive Banach space** for $1 < p < \infty$, the Banach-Alaoglu theorem ensures that any bounded sequence $\{u_k\}$ in $W^{1,p}(\Omega)$ has a subsequence (still denoted by $\{u_k\}$) and a function $u \in W^{1,p}(\Omega)$ such that

$$ u_k \to u \quad \text{in } W^{1,p}(\Omega). \qquad (1.81) $$

This means that for all test functions $\varphi \in W^{1,p'}(\Omega)$, where $p'$ is the Hölder conjugate of $p$ satisfying $\frac{1}{p} + \frac{1}{p'} = 1$, we have

$$ \int_{\Omega} u_k \varphi dx \to \int_{\Omega} u \varphi dx, \quad \int_{\Omega} \nabla u_k \cdot \nabla \varphi dx \to \int_{\Omega} \nabla u \cdot \nabla \varphi dx. \qquad (1.82) $$

However, weak convergence alone does not imply compactness. To obtain strong convergence in $L^q(\Omega)$, we need additional arguments. This is accomplished using the **Fréchet-Kolmogorov**

<figure><img src="image_11.png" /><figcaption>Figure 1.10: **Rellich-Kondrachov Compactness Theorem Illustration**</figcaption></figure>

**compactness criterion**, which states that a bounded subset of $L^q(\Omega)$ is compact if and only if it is **tight** and **uniformly equicontinuous**. More formally, compactness follows if

1. The sequence $u_k(x)$ does not oscillate excessively at small scales.
<figure><img src="image_12.png" /><figcaption>Figure 1.11: L2 norm of $u_n$ converging to 0 (Strong convergence in L2)</figcaption></figure>

2. The sequence $u_k(x)$ does not escape to infinity in a way that prevents strong convergence.

To quantify this, we invoke the **Sobolev-Poincaré inequality**, which states that for $p < n$, there exists a constant $C$ such that

$$ \|u - u_{\Omega}\|_{L^q(\Omega)} \le C \|\nabla u\|_{L^p(\Omega)}, \quad u_{\Omega} = \frac{1}{|\Omega|} \int_{\Omega} u(x) \, dx. \qquad (1.83) $$

Applying this inequality to $u_k - u$, we obtain

$$ \|u_k - u\|_{L^q(\Omega)} \le C \|\nabla(u_k - u)\|_{L^p(\Omega)}. \qquad (1.84) $$

Since $\nabla u_k$ is weakly convergent in $L^p(\Omega)$, we have

$$ \|\nabla u_k - \nabla u\|_{L^p(\Omega)} \to 0. \qquad (1.85) $$

Thus,

$$ \|u_k - u\|_{L^q(\Omega)} \to 0, \qquad (1.86) $$

which establishes the **strong convergence in $L^q(\Omega)$**, completing the proof. The key insight is that compactness arises because the gradients of $u_k$ provide control over the oscillations of $u_k$, ensuring that the sequence cannot oscillate indefinitely without converging in norm. The crucial role of Sobolev embeddings is to guarantee that even though $W^{1,p}(\Omega)$ does not embed compactly into itself, it does embed compactly into $L^q(\Omega)$ for $q < \frac{np}{n-p}$. This embedding ensures that weak convergence in $W^{1,p}(\Omega)$ implies strong convergence in $L^q(\Omega)$, proving the theorem.

### 1.2.4.3 Python Code to generate Figures 1.10 and 1.11 illustrating Rellich-Kondrachov Compactness Theorem

The Python code below produces the Figures 1.10 and 1.11 illustrating Rellich-Kondrachov Compactness Theorem.
```python
import numpy as np
import matplotlib.pyplot as plt

# Define domain
x = np.linspace(0, 1, 500)

# Construct a bounded sequence in H^1(0,1)
# u_n(x) = sin(n*π*x) / n  --> bounded in H^1 but converges strongly in L^2 to 0

def u(n, x):
    return np.sin(n * np.pi * x) / n

# Plot some functions from the sequence
plt.figure(figsize=(10,6))
for n in [1, 2, 5, 10, 20, 50]:
    plt.plot(x, u(n, x), label=f"u_{n}(x)")

plt.title("Rellich-Kondrachov Compactness Theorem Illustration\n"
          "Sequence bounded in H^1 but converging strongly in L^2")
plt.xlabel("x")
plt.ylabel("u_n(x)")
plt.legend()
plt.grid(True)
plt.show()

# Show L2 norm convergence to 0
L2 norms = [np.trapz(u(n, x)**2, x)**0.5 for n in range(1, 51)]

fig, ax = plt.subplots(figsize=(8,5))
ax.tran
ax plot(range(1, 51), L2 norms, marker='o')
ax title("L2 norm of u_n converging to 0 (Strong convergence in L2)")
ax label("n")
ax ylabel(r"$||u_n||_{L^2}$")
ax grid(True)
ax show()
```

## 1.2.5 Fréchet-Kolmogorov compactness criterion

### 1.2.5.1 Literature Review of Fréchet-Kolmogorov compactness criterion

The Fréchet-Kolmogorov compactness criterion is a fundamental result in functional analysis, particularly in the study of compact subsets of $L^p$ spaces, and has been extensively discussed and applied in numerous books and articles. One of the most authoritative references is the classic text by Walter Rudin (1973) [70], which provides a rigorous treatment of the criterion, embedding it within the broader context of weak compactness and equicontinuity in $L^p$ spaces. Rudin's exposition is notable for its clarity and depth, linking the criterion to the Arzelà-Ascoli theorem and emphasizing its utility in proving compactness for families of functions. Another essential reference is by H.L. Royden (1968) [1537], which offers a detailed proof of the Fréchet-Kolmogorov theorem, highlighting its reliance on uniform integrability and tightness conditions. Royden's treatment is particularly valuable for its pedagogical approach, making the criterion accessible to students while maintaining mathematical rigor. The monograph by Robert A. Adams and John J.F. Fournier (2003) [1522] delves into the Fréchet-Kolmogorov criterion in the context of Sobolev embeddings, demonstrating its critical role in establishing compactness for bounded subsets of Sobolev spaces.
Adams and Fournier's work is indispensable for researchers in partial differential equations, as it connects the criterion to applications in variational methods and elliptic regularity theory. Similarly, by Lawrence C. Evans (2010) [1528] discusses the criterion within the framework of weak convergence and compact embeddings, showcasing its importance in the analysis of solutions to PDEs. Evans' presentation is particularly rigorous, with a focus on the interplay between the criterion and the Rellich-Kondrachov theorem. Another key text is by Lawrence C. Evans and Ronald F. Gariepy (1992) [1529], which explores the criterion's implications for functions of bounded variation and its applications in geometric measure theory.

In the realm of probability theory, by Aad W. van der Vaart and Jon A. Wellner (1996) [1540] examines the Fréchet-Kolmogorov criterion as a tool for establishing tightness of probability measures on function spaces. This work is foundational for researchers in empirical process theory, providing a bridge between functional analysis and statistical applications. The authors rigorously develop the criterion's connection to Prokhorov's theorem and its use in proving Donsker-type results. Another probabilistic perspective is offered by Patrick Billingsley (1968) [1524], which discusses the criterion in the context of compactness for stochastic processes. Billingsley's treatment is highly detailed, with a focus on the criterion's role in proving relative compactness for sequences of measures on metric spaces. For a historical perspective, the original works by Maurice Fréchet (1906) [1541] and Andrey Kolmogorov (1931) [1542] are invaluable. Fréchet's early 20th-century papers laid the groundwork for the abstract formulation of compactness in function spaces, while Kolmogorov's contributions refined the criterion, particularly for $L^p$ spaces. These foundational texts are often cited in modern treatments for their pioneering insights. More recently, by Fernando Albiac and Nigel J. Kalton (2006) [1523] provides a modern exposition of the Fréchet-Kolmogorov criterion, embedding it within the broader theory of Banach spaces. Albiac and Kalton's discussion is notable for its emphasis on the criterion's relationship with other compactness principles, such as the Eberlein-Šmulian theorem.

The article by Jürgen Appell and Petr P. Zabrejko (1990) [1543] surveys various compactness criteria, including the Fréchet-Kolmogorov theorem, and discusses their applications in nonlinear analysis. This work is particularly useful for its comprehensive overview and its focus on practical implications for integral operators and fixed-point theorems. Another significant contribution is by Klaus Deimling (1985) [1527], which explores the criterion's use in proving existence theorems for differential equations. Deimling's treatment is rigorous and application-oriented, making it a valuable resource for researchers in nonlinear analysis. In the context of harmonic analysis, by Gerald B. Folland (1992) [1530] discusses the Fréchet-Kolmogorov criterion in relation to compactness for families of oscillatory integrals. Folland's exposition is clear and concise, with a focus on the criterion's role in proving convergence theorems for Fourier series and transforms. Similarly, by Loukas Grafakos (2008) [1531] includes a detailed discussion of the criterion, particularly its applications in the study of singular integrals and maximal functions. Grafakos' work is highly technical but provides deep insights into the criterion's utility in modern harmonic analysis.

The Fréchet-Kolmogorov criterion has also been extensively studied in the context of dynamical systems and ergodic theory. By Karl Petersen (1983) [1536] discusses the criterion's role in proving compactness for families of invariant measures, highlighting its importance in the study of dynamical systems. Petersen's treatment is mathematically rigorous and includes applications to symbolic dynamics and entropy theory. Another relevant reference is by Anatole Katok and Boris Hasselblatt (1995) [1533], which explores the criterion's use in establishing compactness for sequences of diffeomorphisms and its implications for structural stability. In the field of numerical analysis, by K.W. Morton and D.F. Mayers (1994) [1535] discusses the Fréchet-Kolmogorov criterion in the context of finite difference methods. The authors demonstrate how the criterion can be used to prove convergence of numerical schemes by establishing compactness for sequences
of approximate solutions. This practical perspective is complemented by by Philippe G. Ciarlet (1978) [1526], which applies the criterion to prove compactness for finite element approximations. Ciarlet’s work is highly regarded for its rigorous mathematical foundation and its emphasis on error estimation.

The Fréchet-Kolmogorov criterion has also found applications in control theory, as discussed by Jacques-Louis Lions (1971) [1534]. Lions’ work is seminal in this area, showing how the criterion can be used to establish compactness for families of controls and states in PDE-constrained optimization problems. Another important reference is by Roger Temam (1988) [1539], which applies the criterion to prove compactness for attractors in dissipative systems. Temam’s exposition is both deep and broad, covering applications to fluid mechanics and elasticity. Finally, the Fréchet-Kolmogorov criterion has been explored in the context of stochastic PDEs, as seen by Pao-Liu Chow (2007) [1525]. Chow’s work demonstrates the criterion’s utility in proving tightness for solutions to SPDEs, linking it to martingale theory and Skorokhod embeddings. This rigorous treatment is essential for researchers in stochastic analysis. Another notable reference is by Ioannis Karatzas and Steven E. Shreve (1991) [1532], which discusses the criterion’s role in proving compactness for diffusion processes. Karatzas and Shreve’s book is a cornerstone of stochastic analysis, with a thorough and mathematically precise exposition.

### 1.2.5.2 Analysis of Fréchet-Kolmogorov compactness criterion

The Fréchet-Kolmogorov compactness criterion provides necessary and sufficient conditions for a set $\mathcal{F} \subset L^p(\mathbb{R}^n)$, where $1 \le p < \infty$, to be relatively compact, meaning that any sequence $\{f_n\} \subset \mathcal{F}$ has a strongly convergent subsequence in $L^p(\mathbb{R}^n)$. The criterion asserts that $\mathcal{F}$ is relatively compact in $L^p(\mathbb{R}^n)$ if and only if the following three conditions hold: boundedness in $L^p$, uniform integrability (tightness), and equicontinuity in an integral sense. The proof consists of establishing both the necessary and sufficient conditions. To prove necessity, suppose that $\mathcal{F}$ is relatively compact in

<figure><img src="image_13.png" /><figcaption>Figure 1.12: Fréchet–Kolmogorov Criterion: Translation Equicontinuity</figcaption></figure>

$L^p$. Then every sequence $\{f_n\} \subset \mathcal{F}$ has a strongly convergent subsequence. This implies that the
norms of the functions in $\mathcal{F}$ are uniformly bounded, since otherwise there would exist a sequence $\{f_n\}$ such that $\|f_n\|_{L^p} \to \infty$, contradicting compactness. Hence, there exists a constant $M > 0$ such that for all $f \in \mathcal{F}$,

$$ \|f\|_{L^p} = \left( \int_{\mathbb{R}^n} |f(x)|^p dx \right)^{\frac{1}{p}} \leq M. \qquad (1.87) $$

Next, to establish uniform integrability, assume for contradiction that $\mathcal{F}$ is not tight. Then there exists $\varepsilon > 0$ such that for every compact set $K \subset \mathbb{R}^n$, there exists some function $f_K \in \mathcal{F}$ satisfying

$$ \int_{\mathbb{R}^n \setminus K} |f_K(x)|^p dx \geq \varepsilon. \qquad (1.88) $$

This contradicts relative compactness, as it implies the existence of a sequence $\{f_n\}$ with mass escaping to infinity, preventing any strong convergence in $L^p$. Hence, for every $\varepsilon > 0$, there exists a compact $K \subset \mathbb{R}^n$ such that

$$ \sup_{f \in \mathcal{F}} \int_{\mathbb{R}^n \setminus K} |f(x)|^p dx < \varepsilon. \qquad (1.89) $$

To establish the necessity of equicontinuity, suppose for contradiction that it does not hold. Then there exists $\varepsilon > 0$ such that for every $\delta > 0$, there is some function $f_\delta \in \mathcal{F}$ and some shift $h$ with $|h| < \delta$ such that

$$ \int_{\mathbb{R}^n} |f_\delta(x + h) - f_\delta(x)|^p dx \geq \varepsilon. \qquad (1.90) $$

This contradicts compactness, as it implies the existence of arbitrarily oscillatory sequences preventing strong convergence in $L^p$. Hence, for every $\varepsilon > 0$, there exists $\delta > 0$ such that for all $|h| < \delta$,

$$ \sup_{f \in \mathcal{F}} \int_{\mathbb{R}^n} |f(x+h) - f(x)|^p dx < \varepsilon. \qquad (1.91) $$

To prove sufficiency, assume that $\mathcal{F}$ satisfies boundedness, uniform integrability, and equicontinuity. Consider a sequence $\{f_n\} \subset \mathcal{F}$. The first condition guarantees that the sequence is uniformly bounded in $L^p$, ensuring weak compactness by Banach-Alaoglu. The second condition ensures that the functions do not escape to infinity, implying tightness. The third condition ensures that the sequence is equicontinuous, preventing high-frequency oscillations. By the *diagonal argument*, we can extract a subsequence $\{f_{n_k}\}$ that is a Cauchy sequence in $L^p$, ensuring strong convergence. To rigorously show that a subsequence converges, define the *modulus of continuity functional*,

$$ \omega_{\mathcal{F}}(\delta) = \sup_{f \in \mathcal{F}} \int_{\mathbb{R}^n} |f(x+h) - f(x)|^p dx, \qquad (1.92) $$

which satisfies $\omega_{\mathcal{F}}(\delta) \to 0$ as $\delta \to 0$ due to equicontinuity. Given $\varepsilon > 0$, choose $\delta$ such that $\omega_{\mathcal{F}}(\delta) < \frac{\varepsilon}{3}$, and a compact set $K$ such that

$$ \sup_{f \in \mathcal{F}} \int_{\mathbb{R}^n \setminus K} |f(x)|^p dx < \frac{\varepsilon}{3}. \qquad (1.93) $$

By weak compactness, there exists a subsequence $f_{n_k}$ converging weakly in $L^p$ to some $f$. Applying Vitali's theorem, we obtain strong convergence, proving compactness. Thus, the Fréchet-Kolmogorov criterion is fully established.

### 1.2.5.3 Python Code to generate Figure 1.12

The Python code below produces the Figure 1.12 illustrating the Fréchet-Kolmogorov compactness criterion.
```python
import numpy as np
import matplotlib.pyplot as plt

# Domain
x = np.linspace(0, 1, 2000)

# Example 1: Compact candidate family (smooth functions)
def f_family.smooth(n, x):
    return np.sin(np.pi * x) / (1 + n)  # smooth, shrinking amplitude

# Example 2: Non-compact family (high oscillations)
def f_family.oscillatory(n, x):
    return np.sin(n * np.pi * x)  # oscillates faster with n

# L2 norm helper
def L2_norm(f, g, x):
    return npsqrt(np.trapz((f - g)**2, x))

# Translation invariance test
def translationnorms(family, max_n=5, hs= np лinspace(0,0.1,20)):
    norms = []
    for h in hs:
        vals = []
        for n in range(1, max_n+1):
            f = family(n, x)
            f_shift = family(n, (x+h) % 1)  # periodic shift
            vals.append(L2_norm(f, f_shift, x))
        norms.append(max vals))  # sup over family subset
    return hs, norms

# Compute for both families
hs, normsSmooth = translation norms(f_fam ilySmooth)
_, norms_osc = translation norms(f_fam ilyoscid ulatory)

# Plot results
plt. figure(figsize=(8,5))
plt. plot(hs, normsSmooth, marker='o', label="Smooth shrinking family")
plt. plot(hs, norms_osc, marker='s', label="Oscillatory family")
plt. title("Frechet Kolmogorov Criterion Translation Equicontinuity")
plt. xlabel("Shift h")
plt. ylabel("$\sup_n || f(x+h) - f(x) ||_{L^2} $")
plt. legend()
plt. grid(True)
plt. show()
```

## 1.2.6 Sobolev-Poincaré inequality

### 1.2.6.1 Literature Review of Sobolev-Poincaré inequality

The Sobolev-Poincaré inequality is a cornerstone of modern analysis with far-reaching applications across partial differential equations, geometric analysis, and metric geometry. Among the most comprehensive treatments is Vladimir Maz'ya's (1985) [1559] seminal work, which provides an exhaustive study of Sobolev-type inequalities including sharp constants, capacity theory, and their
validity on both smooth and irregular domains. This foundational text is complemented by Haim Brezis's (2011) [1549] contribution, which situates the inequality within functional analysis and PDE theory, particularly emphasizing its role in compact embeddings and elliptic regularity. The geometric aspects are profoundly explored in Mikhail Gromov's (1999) [1555] work, extending classical results to general metric measure spaces, while Elliott Lieb and Michael Loss's (2001) [1558] analysis offers crucial insights into optimal constants and connections with isoperimetric inequalities through real analytic methods. Heinonen's (2001)[1563] lectures on analysis on metric spaces also discuss the Sobolev-Poincaré inequality and its application in various mathematical topics in great depth.

For applications to partial differential equations, David Gilbarg and Neil Trudinger's (2001) [1562] treatment remains indispensable for its rigorous examination of local and global versions on bounded domains, particularly in regularity theory. The variational perspective appears prominently in Antonio Ambrosetti and Andrea Malchiodi's (2007) [1546] study, which demonstrates how the inequality underpins existence results for nonlinear PDEs. In geometric contexts, Pierre Buser's (1982) [1550] research on Riemannian manifolds establishes precise relationships between Poincaré constants and geometric data like curvature and diameter, whereas Nicola Garofalo and Dariusz Z. D'Agostino's (1999) [1554] contributions derive sharp inequalities for sub-Riemannian structures. The probabilistic and geometric-analytic connections are masterfully presented in Dominique Bakry, Ivan Gentil, and Michel Ledoux's (2014) [1547] work, linking the inequality to Markov semigroups and curvature-dimension conditions.

The extension to nonsmooth settings is systematically developed in Juha Heinonen and Pekka Koskela's (2001, 2005) [1556][1557] research, which reformulates the inequality using upper gradients for applications in fractal and sub-Riemannian geometry. Luigi Ambrosio, Nicola Gigli, and Giuseppe Savaré's (2005, 2008) [1544][1545] contributions provide groundbreaking syntheses with optimal transport and metric measure space theory. Hypoelliptic settings are addressed in Luca Capogna, Donatella Danielli, and Nicola Garofalo's (2001) [1551] work on Hörmander vector fields, while Fabrice Baudoin and Michel Bonnefont's (2012) [1548] research extends these to subelliptic operators. The heat kernel and spectral perspectives are thoroughly examined in E. Brian Davies's (1989) [1552] investigation and Qi S. Zhang's (2011) [1561] work on Ricci flow. Further abstract formulations appear in Karl-Theodor Sturm and Max-K. von Renesse's (2005) [1560] optimal transport approach to curvature bounds, and Herbert Federer's (1969) [1553] study provides essential links to currents and varifolds. Collectively, these works form a multifaceted theoretical edifice, demonstrating the Sobolev-Poincaré inequality's central role in analysis, geometry, and probability through diverse yet interconnected perspectives.

### 1.2.6.2 Analysis of Sobolev-Poincaré inequality

Let $\Omega \subset \mathbb{R}^n$ be a bounded domain with a Lipschitz boundary. Consider the Sobolev space $W^{1,p}(\Omega)$ for $1 \le p < n$, which consists of functions $u \in L^p(\Omega)$ whose weak derivatives $\nabla u$ also belong to $L^p(\Omega)$. The Sobolev-Poincaré inequality asserts the existence of a constant $C = C(\Omega, p) > 0$ such that

$$ \|u - u_{\Omega}\|_{L^{p^{*}}(\Omega)} \le C \|\nabla u\|_{L^{p}(\Omega)} \quad (1.94) $$

for all $u \in W^{1,p}(\Omega)$, where $p^* = \frac{np}{n-p}$ is the Sobolev conjugate exponent and $u_{\Omega}$ is the mean value of $u$ over $\Omega$, given by

$$ u_{\Omega} = \frac{1}{|\Omega|} \int_{\Omega} u(x) dx. \quad (1.95) $$

To prove this inequality, we first use a representation formula for $u(x)$. Given any two points $x, y \in \Omega$, we consider the fundamental theorem of calculus applied along the segment connecting $x$
to $y$. Define the parametrization

$$ \gamma(t) = x + t(y - x), \quad t \in [0, 1]. \tag{1.96} $$

Applying the fundamental theorem of calculus along this line segment,

$$ u(y) - u(x) = \int_{0}^{1} \frac{d}{dt} u(\gamma(t)) dt. \tag{1.97} $$

By the chain rule,

$$ \frac{d}{dt} u(\gamma(t)) = \nabla u(\gamma(t)) \cdot (y - x). \tag{1.98} $$

Substituting this into the integral expression,

$$ u(y) - u(x) = \int_{0}^{1} \nabla u(\gamma(t)) \cdot (y - x) dt. \tag{1.99} $$

Taking the absolute value and applying Hölder's inequality with conjugate exponents $p$ and $p' = \frac{p}{p-1}$,

$$ |u(y) - u(x)| \le |y - x| \int_{0}^{1} |\nabla u(\gamma(t))| dt. \tag{1.100} $$

Now, integrating both sides over $y \in \Omega$, using Fubini's theorem and Minkowski's integral inequality,

$$ \int_{\Omega} |u(y) - u(x)|^{p^{*}} dy \le \int_{\Omega} \left( |y - x| \int_{0}^{1} |\nabla u(\gamma(t))| dt \right)^{p^{*}} dy. \tag{1.101} $$

Using the properties of integral norms and applying Hölder's inequality again,

$$ \left( \int_{\Omega} |u(y) - u(x)|^{p^{*}} dy \right)^{\frac{1}{p^{*}}} \le C \|\nabla u\|_{L^p(\Omega)}. \tag{1.102} $$

By averaging over $x \in \Omega$ and defining the mean-value term appropriately, we obtain

$$ \|u - u_{\Omega}\|_{L^{p^{*}}(\Omega)} \le C \|\nabla u\|_{L^{p}(\Omega)}, \tag{1.103} $$

where $C$ depends on $\Omega$ and $p$, completing the proof. Let $\Omega \subset \mathbb{R}^n$ be a bounded domain with a Lipschitz boundary. Consider the Sobolev space $W^{1,p}(\Omega)$ for $1 \le p < n$, which consists of functions $u \in L^p(\Omega)$ whose weak derivatives $\nabla u$ also belong to $L^p(\Omega)$. The Sobolev-Poincaré inequality asserts the existence of a constant $C = C(\Omega, p) > 0$ such that

$$ \|u - u_{\Omega}\|_{L^{p^{*}}(\Omega)} \le C \|\nabla u\|_{L^{p}(\Omega)} \tag{1.104} $$

for all $u \in W^{1,p}(\Omega)$, where $p^* = \frac{np}{n-p}$ is the Sobolev conjugate exponent and $u_{\Omega}$ is the mean value of $u$ over $\Omega$, given by

$$ u_{\Omega} = \frac{1}{|\Omega|} \int_{\Omega} u(x) dx. \tag{1.105} $$

To prove this inequality, we first use a representation formula for $u(x)$. Given any two points $x, y \in \Omega$, we consider the fundamental theorem of calculus applied along the segment connecting $x$ to $y$. Define the parametrization

$$ \gamma(t) = x + t(y - x), \quad t \in [0, 1]. \tag{1.106} $$

Applying the fundamental theorem of calculus along this line segment,

$$ u(y) - u(x) = \int_{0}^{1} \frac{d}{dt} u(\gamma(t)) dt. \tag{1.107} $$
<figure><img src="image_14.png" /><figcaption>Figure 1.13: Test Functions in $H_0^1(0, 1)$</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 1.14: Sobolev-Poincaré Inequality Ratios: $||u||_{L^2} / ||u'||_{L^2}$</figcaption></figure>

By the chain rule,

$$ \frac{d}{dt}u(\gamma(t)) = \nabla u(\gamma(t)) \cdot (y - x). \qquad (1.108) $$

Substituting this into the integral expression,

$$ u(y) - u(x) = \int_{0}^{1} \nabla u(\gamma(t)) \cdot (y - x) dt. \qquad (1.109) $$
Taking the absolute value and applying Hölder's inequality with conjugate exponents $p$ and $p' = \frac{p}{p-1}$,

$$|u(y) - u(x)| \le |y - x| \int_0^1 |\nabla u(\gamma(t))| dt. \qquad (1.110)$$

Now, integrating both sides over $y \in \Omega$, using Fubini's theorem and Minkowski's integral inequality,

$$\int_{\Omega} |u(y) - u(x)|^{p^*} dy \le \int_{\Omega} \left( |y - x| \int_{0}^{1} |\nabla u(\gamma(t))| dt \right)^{p^*} dy. \qquad (1.111)$$

Using the properties of integral norms and applying Hölder's inequality again,

$$\left( \int_{\Omega} |u(y) - u(x)|^{p^*} dy \right)^{\frac{1}{p^*}} \le C \|\nabla u\|_{L^p(\Omega)}. \qquad (1.112)$$

By averaging over $x \in \Omega$ and defining the mean-value term appropriately, we obtain

$$\|u - u_{\Omega}\|_{L^{p^{*}}(\Omega)} \le C \|\nabla u\|_{L^{p}(\Omega)}, \qquad (1.113)$$

where $C$ depends on $\Omega$ and $p$, completing the proof. The consequences and applications are:

1. **Regularity of PDE Solutions:** The Sobolev-Poincaré inequality is crucial in proving the existence and regularity of weak solutions to elliptic PDEs.

2. **Compactness and Rellich-Kondrachov Theorem:** It plays a role in proving the compact embedding of $W^{1,p}(\Omega)$ into $L^q(\Omega)$, which is fundamental in functional analysis.

3. **Control of Function Oscillations:** It quantifies how much a function can deviate from its mean, which is used in various areas of mathematical physics and geometry.

One important extension is the sharp form of the Sobolev-Poincaré inequality, which involves explicit best constants in terms of the domain geometry. Specifically, in some cases, the optimal constant is related to eigenvalues of the Laplace operator or geometric properties such as the diameter or inradius of $\Omega$. Another important extension is the fractional Sobolev-Poincaré inequality, which deals with function spaces like $W^{s,p}(\Omega)$ where $s$ is fractional, incorporating nonlocal effects.

In conclusion, the Sobolev-Poincaré inequality is a powerful mathematical result that provides a precise relationship between function averages and their gradients. It is a cornerstone in modern analysis, with deep implications in PDE theory, functional analysis, and geometric measure theory.

### 1.2.6.3 Python Code to generate Figures 1.13 and 1.14

The Python code below produces the Figures 1.13 and 1.14 illustrating the Fréchet-Kolmogorov compactness criterion.

```python
import numpy as np
import matplotlib.pyplot as plt

# Domain
def grid(N=2000):
    return np.linspace(0, 1, N)

x = grid()

# Define some H_0^1(0,1) test functions (vanish at boundary)
def u1(x): return np.sin(np.pi * x)
```
```python
def u2(x): return x * (1-x)
def u3(x): return np.sin(2*np.pi*x) * (1-x)
def u4(x): return np Paght(x*(1-x)) * np.sine(np.pi*x)

functions = [u1, u2, u3, u4]
labels = ["sin(π x)", "x(1-x)", "sin(2π x)(1-x)", "sqrt(x(1-x)) sin(π x)"]

# L2 norm helper
def L2_norm(f, x):
    return np .sqrt(np.trapz(f(x)**2, x))

# Derivative approximation
def derivative(f, x):
    h = x[1] - x[0]
    return np .gradient(f(x), h)

ratios = []
for f in functions:
    norm_u = L2_norm(f, x)
    norm_du = L2_norm(lambda z: derivative(f, x), x)
    ratios.append norm_u / norm_du)

# Plot functions
plt .figure(figure size=(10,6))
for f, label in zip(functions, labels):
    plt .plot(x, f(x), label =label)
    plt .title("Test Functions in H_0^1(0,1)")
    plt .xlabel("x")
    plt .ylabel("u(x)")
    plt .legend()
    plt .grid(True)
    plt .show()

# Plot ratios ||u|| / ||u'||
45 plt .figure(figure size=(8,5))
46 plt .bar labels, ratios
47
48 plt .title("Sobolev Poincare Inequality Ratios")
49
49 plt .ylabel("Ratio")
50
51 plt .grid(True, axis = 'y')
51 plt .show()
```
## 2 Universal Approximation Theorem: Re-fined Proof

The Universal Approximation Theorem (UAT) is a fundamental result in neural network theory, stating that a feedforward neural network with a single hidden layer containing a finite number of neurons can approximate any continuous function on a compact subset of $\mathbb{R}^n$ to any desired degree of accuracy, provided that an appropriate activation function is used. This theorem has significant implications in machine learning, function approximation, and deep learning architectures.

The Universal Approximation Theorem (UAT) is a fundamental result in the theory of artificial neural networks, which states that a feedforward neural network with a single hidden layer and a non-polynomial activation function can approximate any continuous function on a compact subset of $\mathbb{R}^n$ to arbitrary precision, provided the hidden layer has a sufficient number of neurons. Formally, let $\sigma: \mathbb{R} \to \mathbb{R}$ be a non-constant, bounded, and continuous activation function. For any continuous function $f: K \subset \mathbb{R}^n \to \mathbb{R}$ defined on a compact set $K$, and for every $\epsilon > 0$, there exists a neural network $F: \mathbb{R}^n \to \mathbb{R}$ of the form:

$$F(\mathbf{x}) = \sum_{i=1}^{N} v_i \sigma(\mathbf{w}_i^T \mathbf{x} + b_i), \qquad (2.1)$$

where $N$ is the number of hidden units, $\mathbf{w}_i \in \mathbb{R}^n$ are the input weights, $b_i \in \mathbb{R}$ are the biases, and $v_i \in \mathbb{R}$ are the output weights, such that:

$$\sup_{\mathbf{x} \in K} |F(\mathbf{x}) - f(\mathbf{x})| < \epsilon. \qquad (2.2)$$

The theorem underscores the expressive power of neural networks, demonstrating that even a single hidden layer is theoretically capable of approximating any continuous function, given enough neurons. The key assumptions are the non-polynomial nature of the activation function and the compactness of the domain. Common activation functions satisfying these conditions include the sigmoid $\sigma(x) = \frac{1}{1+e^{-x}}$, the hyperbolic tangent $\tanh(x)$, and the rectified linear unit (ReLU) $\max(0, x)$. The compactness requirement ensures that the function $f$ is uniformly continuous, enabling the approximation to hold uniformly across the entire domain. The proof of the UAT typically relies on tools from functional analysis, such as the Stone-Weierstrass theorem or the Hahn-Banach theorem, which guarantee the density of certain function spaces. For instance, the Stone-Weierstrass theorem asserts that any subalgebra of continuous functions that separates points and contains the constant functions is dense in the space of continuous functions. By showing that the set of functions representable by a neural network forms such a subalgebra, the UAT can be derived. The non-polynomial condition on $\sigma$ is critical to ensure the network can generate a rich enough class of functions to approximate arbitrary continuous mappings.

However, the theorem does not provide explicit bounds on the number of neurons $N$ required to achieve a given accuracy $\epsilon$, nor does it address the practical challenges of training such networks. In practice, deeper architectures with multiple hidden layers often outperform shallow networks, even though the UAT theoretically guarantees the sufficiency of a single hidden layer. This discrepancy
highlights the difference between existence and learnability: while the UAT ensures the existence of an approximating network, finding such a network through optimization remains a non-trivial task. The theorem also does not account for the curse of dimensionality, which can make the required number of neurons infeasibly large for high-dimensional inputs. Despite these limitations, the UAT remains a cornerstone of neural network theory, providing a rigorous foundation for their universal applicability.

## 2.1 Literature Review of Universal Approximation Theorem

Hornik et. al. (1989) [58] in their seminal paper rigorously proved that multilayer feedforward neural networks with a single hidden layer and a sigmoid activation function can approximate any continuous function on a compact set. It extends prior results and lays the foundation for the modern understanding of UAT. Cybenko (1989) [60] provided one of the first rigorous proofs of the UAT using the sigmoid function as the activation function. They demonstrated that a single hidden layer network can approximate any continuous function arbitrarily well. Barron (1993) [61] extended UAT by quantifying the approximation error and analyzing the rate of convergence. This work is crucial for understanding the practical efficiency of neural networks. Pinkus (1999) [62] provided a comprehensive survey of UAT from the perspective of approximation theory and also discussed conditions for approximation with different activation functions and the theoretical limits of neural networks. Lu et.al. (2017) [63] investigated how the width of neural networks affects their approximation capability, challenging the notion that deeper networks are always better. They also provided insights into trade-offs between depth and width. Hanin and Sellke (2018) [64] extended UAT to ReLU activation functions, showing that deep ReLU networks achieve universal approximation while maintaining minimal width constraints. Garcia-Cervera et. al. (2024) [66] extended the universal approximation theorem to set-valued functions and its applications to Deep Operator Networks (DeepONets), which are useful in control theory and PDE modeling. Majee et.al. (2024) [67] explored the universal approximation properties of deep neural networks for solving inverse problems using Markov Chain Monte Carlo (MCMC) techniques. Toscano et. al. (2024) [68] introduced Kurkova-Kolmogorov-Arnold Networks (KKANs), an extension of UAT incorporating Kolmogorov's superposition theorem for improved approximation capabilities. Son (2025) [69] established a new framework for operator learning based on the UAT, providing a theoretical foundation for backpropagation-free deep networks.

## 2.2 Approximation Using Convolution Operators

Let us begin by considering the convolution operator and its role in approximating functions in the context of the Universal Approximation Theorem (UAT). Suppose $f: \mathbb{R}^n \to \mathbb{R}$ is a continuous and bounded function. The convolution of $f$ with a kernel function $\phi: \mathbb{R}^n \to \mathbb{R}$, denoted as $f * \phi$, is defined as

$$ (f * \phi)(x) = \int_{\mathbb{R}^n} f(y)\phi(x - y) dy. \qquad (2.3) $$

The kernel $\phi(x)$ is typically chosen to be smooth, compactly supported, and normalized such that

$$ \int_{\mathbb{R}^n} \phi(x) dx = 1. \qquad (2.4) $$

To approximate $f$ locally, we introduce a scaling parameter $\epsilon > 0$ and define the scaled kernel $\phi_{\epsilon}(x)$ as

$$ \phi_{\epsilon}(x) = \epsilon^{-n} \phi\left(\frac{x}{\epsilon}\right). \qquad (2.5) $$
The factor $\epsilon^{-n}$ ensures that $\phi_{\epsilon}(x)$ remains a probability density function, satisfying

$$ \int_{\mathbb{R}^n} \phi_{\epsilon}(x) dx = \int_{\mathbb{R}^n} \phi(x) dx = 1. \qquad (2.6) $$

The convolution of $f$ with the scaled kernel $\phi_{\epsilon}$ is given by

$$ (f * \phi_{\epsilon})(x) = \int_{\mathbb{R}^n} f(y) \phi_{\epsilon}(x - y) dy. \qquad (2.7) $$

Performing the change of variables $z = \frac{x-y}{\epsilon}$, we have $y = x - \epsilon z$ and $dy = \epsilon^n dz$. Substituting into the integral, we obtain

$$ (f * \phi_{\epsilon})(x) = \int_{\mathbb{R}^n} f(x - \epsilon z) \phi(z) dz. \qquad (2.8) $$

This representation shows that $(f * \phi_{\epsilon})(x)$ is a smoothed version of $f(x)$, where the smoothing is controlled by the parameter $\epsilon$. As $\epsilon \to 0$, the kernel $\phi_{\epsilon}(x)$ becomes increasingly concentrated around $x$, and we recover $f(x)$ in the limit:

$$ \lim_{\epsilon \to 0} (f * \phi_{\epsilon})(x) = f(x), \qquad (2.9) $$

assuming $f$ is continuous. This result can be rigorously proven using properties of the kernel

<figure><img src="image_16.png" /><figcaption>Figure 2.1: Universal Approximation Theorem Illustration</figcaption></figure>

$\phi_{\epsilon}$, such as its smoothness and compact support, and the dominated convergence theorem, which ensures that the integral converges uniformly to $f(x)$. Now, let us consider the role of convolution
operators in the approximation of $f$ by neural networks. A single-layer feedforward neural network is expressed as

$$ \hat{f}(x) = \sum_{i=1}^{M} c_i \sigma(w_i^T x + b_i), \qquad (2.10) $$

where $c_i \in \mathbb{R}$ are coefficients, $w_i \in \mathbb{R}^n$ are weight vectors, $b_i \in \mathbb{R}$ are biases, and $\sigma: \mathbb{R} \to \mathbb{R}$ is the activation function. The activation function $\sigma(w_i^T x + b_i)$ can be interpreted as a localized response function, analogous to the kernel $\phi(x - y)$ in convolution. By drawing an analogy between the two, we can write the neural network approximation as

$$ \hat{f}(x) \approx \sum_{i=1}^{M} f(x_i) \phi_{\epsilon}(x - x_i) \Delta x \qquad (2.11) $$

where $\phi_{\epsilon}(x)$ is interpreted as a parameterized kernel defined by $w_i$, $b_i$, and $\sigma$, and $\Delta x$ represents a discretization step. The approximation error $\|f - \hat{f}\|_{\infty}$ can be decomposed into two components:

$$ \|f - \hat{f}\|_{\infty} \le \|f - f * \phi_{\epsilon}\|_{\infty} + \|f * \phi_{\epsilon} - \hat{f}\|_{\infty}. \qquad (2.12) $$

The term $\|f - f * \phi_{\epsilon}\|_{\infty}$ represents the error introduced by smoothing $f$ with the kernel $\phi_{\epsilon}$, and it can be made arbitrarily small by choosing $\epsilon$ sufficiently small, provided $f$ is regular enough (e.g., Lipschitz continuous). The term $\|f * \phi_{\epsilon} - \hat{f}\|_{\infty}$ quantifies the error due to discretization, which

<figure><img src="image_17.png" /><figcaption>Figure 2.2: Approximation Error vs Hidden Layer Size</figcaption></figure>

vanishes as the number of neurons $M \to \infty$. To rigorously analyze the convergence of $\hat{f}(x)$ to $f(x)$, we rely on the density of neural network approximators in function spaces. The Universal
Approximation Theorem states that, for any continuous function $f$ on a compact domain $\Omega \subset \mathbb{R}^n$ and any $\epsilon > 0$, there exists a neural network $\hat{f}$ with finitely many neurons such that

$$ \sup_{x \in \Omega} |f(x) - \hat{f}(x)| < \epsilon. \qquad (2.13) $$

This result hinges on the ability of the activation function $\sigma$ to generate a rich set of basis functions. For example, if $\sigma(x) = \max(0, x)$ (ReLU), the network approximates $f(x)$ by piecewise linear functions. If $\sigma(x) = \frac{1}{1+e^{-x}}$ (sigmoid), the network generates smooth approximations that resemble logistic regression.

In this refined proof of the UAT, convolution operators provide a unifying framework for understanding the smoothing, localization, and discretization processes that underlie neural network approximations. The interplay between $\phi_{\epsilon}(x)$, $f * \phi_{\epsilon}(x)$, and $\hat{f}(x)$ reveals the profound mathematical structure that connects classical approximation theory with modern machine learning. This connection not only enhances our theoretical understanding of neural networks but also guides the design of architectures and algorithms for practical applications.

### 2.2.1 Python Code to generate Figures 2.1 and 2.2

The Python code below produces the Figures 2.1 and 2.2 illustrating the Approximation Using Convolution Operators.

```python
import torch
import torch.nn as nn
import torch.optim as optim
import matplotlib.pyplot as plt
import numpy as np

# Target function
def f(x):
    return np.sin(2 * np.pi * x)

# Data
x np = np.linspace(0, 1, 200)
y np = f(x np)

x = torch.tensor(x np, dtype = torch.float32).view(-1, 1)
y = torch.tensor(y np, dtpe = torch.float32).view(-1, 1)

# Define simple feedforward network with variable hidden size
class Net (nn.Module):
    def __init__ (self, hidden_size):
        super (Net, self).__init__ ()
        self hidden = nn.Linear(1, hidden_size)
        self activation = nn.Tanh() # non-linear activation
        self output = nn.Linear(hidden_size, 1)

    def forward(self, x):
        return self .output(self .activation(self .hidden(x)))

# Train and plot approximations for different hidden sizes
hidden_size s = [2, 5, 10, 20, 50]
errors = []

plt .figure (figsize=(12, 8))
```
```python
34
35 for i, h in enumerate(hidden_sizes, 1):
36     net = Net(h)
37     criterion = nn.MSELoss()
38     optimizer = optim.Adam(net.parameters(), lr=0.05)
39
40 # Train
41 for epoch in range(2000):
42     optimizer.zeroGrad()
43     y_pred = net(x)
44     loss = criterion(y_pred, y)
45     loss backward()
46     optimizer.step()
47
48 # Evaluate
49 with torch.noGrad():
50     y_fit = net(x).numpy()
51     error = npsqrt(np.trapz((y_fit.flatten() - y NP)**2, x NP))
52     errors.append error)
53
54 plt Subplot(2, 3, i)
55 plt. plot(x np, y np, 'b-', label="Target sin(2πx)")
56 plt. plot(x np, y fit, 'r--', label="NN approx (hidden={h})")
57 plt. title(f"Hidden layer size = {h}")
58 plt. legend()
59 plt. grid(True)
60
61 plt. subTitle("Universal Approximation Theorem Illustration\n"
62               "Single hidden-layer NN approximating sin(2πx)")
63 plt. tight layout (rect=[0, 0, 1, 0.96])
64 plt. show()
65
66 # Plot error vs hidden size
67 plt. figure (figsize=(7,5))
68 plt. plot(hidden sizes, errors, marker='o')
69 plt. title("Approximation Error vs Hidden Layer Size")
70 plt.xlabel("Hidden layer size")
71 plt. ylabel("$L^2$ error")
72 plt. grid (True)
73 plt. show()
```

### 2.2.2 Stone-Weierstrass Application

#### 2.2.2.1 Literature Review of Stone-Weierstrass Application

Rudin (1976) [70] introduced the Weierstrass approximation theorem and proves its generalization, the Stone-Weierstrass theorem. He also discussed the algebraic structure of function spaces and how the theorem ensures the uniform approximation of continuous functions by polynomials. He also presented examples and exercises related to compactness, uniform convergence, and Banach algebra structures. Stein and Shakarchi (2005) [71] extended the Stone-Weierstrass theorem into measure theory and functional analysis. He also proved the theorem in the context of Lebesgue integration. He also discussed how it applies to Hilbert spaces and orthogonal polynomials. He also connected the theorem to Fourier analysis and spectral decomposition. Conway (2019) [73] explored the Stone-Weierstrass theorem in the setting of Banach algebras and C-algebras*. He also extended
the theorem to non-commutative function algebras and discussed the operator-theoretic implications of the theorem in Hilbert spaces. He also analyzed the theorem's application to spectral theory. Dieudonné (1981) [74] traced the historical development of functional analysis, including the origins of the Stone-Weierstrass theorem and discussed contributions by Karl Weierstrass and Marshall Stone. He also explored how the theorem influenced topological vector spaces and operator theory and also included perspectives on the axiomatic development of function approximation. Folland (1999) [76] discussed the Stone-Weierstrass theorem in depth with applications to probability theory and ergodic theory and used the theorem to establish the density of algebraic functions in measure spaces. He also connected the Stone-Weierstrass theorem to functional approximation in Lp spaces. He also explored the interplay between the Stone-Weierstrass theorem and the Hahn-Banach theorem. Sugiura (2024) [77] extended the Stone-Weierstrass theorem to the study of reservoir computing in machine learning and proved that certain neural networks can approximate functions uniformly under the assumptions of the theorem. He bridges classical functional approximation with modern AI and deep learning. Liu et al. (2024) [78] investigated the Stone-Weierstrass theorem in normed module settings and used category theory to generalize function approximation results. He also extended the theorem beyond real-valued functions to structured mathematical objects. Martinez-Barreto (2025) [79] provided a modern formulation of the theorem with rigorous proof and reviewed applications in operator algebras and topology. He also discussed open problems related to function approximation. Chang and Wei (2024) [80] used the Stone-Weierstrass theorem to derive new operator inequalities and applied the theorem to functional analysis in quantum mechanics. Caballer et al. (2024) [81] investigated cases where the Stone-Weierstrass theorem fails and provided counterexamples and refined conditions for uniform approximation. Chen (2024) [82] extended the Stone-Weierstrass theorem to generalized function spaces and introduced a new class of uniform topological algebras. Rafiei and Akbarzadeh-T (2024) [83] used the Stone-Weierstrass theorem to analyze function approximation in fuzzy logic systems and explored the applications in control systems and AI.

### 2.2.2.2 Analysis of Stone-Weierstrass Application

The **Stone-Weierstrass Theorem** serves as a cornerstone in functional analysis, bridging the algebraic structure of continuous functions with approximation theory. This theorem, when applied to the **Universal Approximation Theorem (UAT)**, provides a rigorous foundation for asserting that neural networks can approximate any continuous function defined on a compact set. To understand this connection in its most scientifically and mathematically rigorous form, we must carefully analyze the algebra of continuous functions on a compact Hausdorff space and the role of neural networks in approximating these functions, ensuring that all mathematical nuances are explored with extreme precision. Let $X$ be a compact Hausdorff space, and let $C(X)$ represent the space of continuous real-valued functions on $X$. The **supremum norm** $\|f\|_{\infty}$ for a function $f \in C(X)$ is defined as:

$$ \|f\|_{\infty} = \sup_{x \in X} |f(x)| \qquad (2.14) $$

This supremum norm is critical in defining the proximity between continuous functions, as we seek to approximate any function $f \in C(X)$ by a function $g$ from a subalgebra $A \subset C(X)$. The **Stone-Weierstrass theorem** guarantees that if the subalgebra $A$ satisfies two essential properties—(1) it contains the constant functions, and (2) it separates points—then the closure of $A$ in the supremum norm will be the entire space $C(X)$. To formalize this, we define the **point separation property** as follows: for every pair of distinct points $x_1, x_2 \in X$, there exists a function $h \in A$ such that $h(x_1) \neq h(x_2)$. This condition ensures that functions from $A$ are sufficiently “rich” to distinguish between different points in $X$. Mathematically, this is expressed as:

$$ \exists h \in A \text{ such that } h(x_1) \neq h(x_2) \quad \forall x_1, x_2 \in X, x_1 \neq x_2 \qquad (2.15) $$
<figure><img src="image_18.png" /><figcaption>Figure 2.3: Stone–Weierstrass (Weierstrass) demo: Bernstein polynomial approximations</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 2.4: Zoom near the kink ($x = 0.3$) — convergence is uniform</figcaption></figure>

Given these two properties, the Stone-Weierstrass theorem asserts that for any continuous function $f \in C(X)$ and any $\epsilon > 0$, there exists an element $g \in A$ such that:

$$ \|f - g\|_{\infty} < \epsilon \qquad (2.16) $$

This result ensures that any continuous function on a compact Hausdorff space can be approximated arbitrarily closely by functions from a sufficiently rich subalgebra. In the context of the **Universal Approximation Theorem (UAT)**, we seek to apply the Stone-Weierstrass theorem to the approximation capabilities of neural networks. Let $K \subseteq \mathbb{R}^n$ be a compact subset, and let $f \in C(K)$ be a continuous function defined on this set. A feedforward neural network with a non-linear activation function $\sigma$ has the form:

$$ \hat{f}_{\theta}(x) = \sum_{i=1}^{N} w_i \sigma(\langle \mathbf{w}_i, x \rangle + b_i) \qquad (2.17) $$
where $\langle \mathbf{w}_i, x \rangle$ represents the inner product between the weight vector $\mathbf{w}_i$ and the input $x$, and $b_i$ represents the bias term. The activation function $\sigma$ is typically non-linear (such as the sigmoid or ReLU function), and the parameters $\theta = \{w_i, b_i\}_{i=1}^N$ are the weights and biases of the network. The function $\hat{f}_{\theta}(x)$ is a weighted sum of the non-linear activations applied to the affine transformations of $x$.

We now explore the connection between neural networks and the Stone-Weierstrass theorem. A critical observation is that the set of functions defined by a neural network with non-linear activation is a subalgebra of $C(K)$ provided the activation function $\sigma$ is sufficiently rich in its non-linearity. This non-linearity ensures that the network can separate points in $K$, meaning that for any two

<figure><img src="image_20.png" /><figcaption>Figure 2.5: **Sup-norm error** $||f - B_n(f)||_{\infty}$ (log-log scale)</figcaption></figure>

distinct points $x_1, x_2 \in K$, there exists a network function $\hat{f}_{\theta}$ that takes distinct values at these points. This satisfies the point separation condition required by the Stone-Weierstrass theorem. To formalize this, consider two distinct points $x_1, x_2 \in K$. Since $\sigma$ is non-linear, the function $\hat{f}_{\theta}(x)$ with appropriately chosen weights and biases will satisfy:

$$ \hat{f}_{\theta}(x_1) \neq \hat{f}_{\theta}(x_2) \qquad (2.18) $$

Thus, the algebra of neural network functions satisfies the point separation property. By applying the Stone-Weierstrass theorem, we conclude that this algebra is dense in $C(K)$, meaning that for any continuous function $f \in C(K)$ and any $\epsilon > 0$, there exists a neural network function $\hat{f}_{\theta}$ such that:

$$ \|f(x) - \hat{f}_{\theta}(x)\|_{\infty} < \epsilon \quad \forall x \in K \qquad (2.19) $$

This rigorous result shows that neural networks with a non-linear activation function can approximate any continuous function on a compact set arbitrarily closely in the supremum norm, thereby proving the Universal Approximation Theorem. To further explore this, consider the error term:

$$ \|f(x) - \hat{f}_{\theta}(x)\|_{\infty} \qquad (2.20) $$
For a given function *f* and a compact set *K*, this error term can be made arbitrarily small by increasing the number of neurons in the hidden layer of the neural network. This increases the capacity of the network, effectively enlarging the subalgebra of functions generated by the network, thereby improving the approximation. As the number of neurons increases, the network's ability to approximate any function from *C*(*K*) becomes increasingly precise, which aligns with the conclusion of the Stone-Weierstrass theorem that the network functions form a dense subalgebra in *C*(*K*). Thus, the Universal Approximation Theorem, derived through the Stone-Weierstrass theorem, rigorously proves that neural networks can approximate any continuous function on a compact set to any desired degree of accuracy. The combination of the non-linearity of the activation function and the architecture of the neural network guarantees that the network can generate a dense subalgebra of continuous functions, ultimately allowing it to approximate any function from *C*(*K*). This result not only formalizes the approximation power of neural networks but also provides a deep theoretical foundation for understanding their capabilities as universal approximators.

### 2.2.2.3 Python Code to generate Figures 2.3, 2.4 and 2.5

The Python code below produces the Figures 2.3, 2.4 and 2.5 illustrating the Stone-Weierstrass Application.

```python
"""
Visual demonstration of the Stone-Weierstrass (Weierstrass) theorem via
Bernstein polynomial approximation on [0,1].
"""

import numpy as np
import matplotlib.pyplot as plt
import math
from scipy.special import comb  # for binomial coefficients

# ---------- Choose a continuous target function on [0,1] ----------
def f(x):
    # A continuous but non-differentiable example (kink at x=0.3)
    return np.abs(x - 0.3)

# alternative examples:
def f(x): return np.sin(2*np.pi*x) + 0.5*x**2
# def f(x): return npwhere(x<0.5, x, 1-x) # continuous triangle-like

# dense grid for plotting and sup-norm estimation
X = np.linspace(0.0, 1.0, 2001)
fX = f(X)

# ---------- Bernstein polynomial B_n(f) ----------
def bernsteinApprox(f vals_at_nodes, n, x_points):
    """
    Compute Bernstein polynomial of degree n using values f(k/n)
    provided in f vals_at_nodes.
    f vals_at_nodes should be array of length n+1: [f(0), f(1/n), ..., f
    (1)].
    x_points may be a numpy array.
    Returns B_n(f)(x_points).
    """
    x = x_points
    # Use vectorized evaluation: sum_{k=0}^n f(k/n) * C(n,k) x^k (1-x)^{
```
```python
k = np.arange(n+1)
# shape (n+1, len(x))
# compute basis B_{n,k}(x) for each k
# to avoid overflow for large n we compute in logspace if needed,
# but for moderate n ( <= 500) direct comb is fine.
basis = np.array([comb(n, ki) * (x**ki) * ((1-x)**(n-ki)) for ki in
                  k])
# weighted sum
return np.dot(f vals at nodes, basis)
```

```python
# ---------- Compute approximations for several n ----------
degrees = [2, 5, 10, 20, 50, 100] # πck some n's to show
approxs = {}
sup_errors = []
```

```python
for n in degrees:
    nodes = np.:linspace(0.0, 1.0, n+1)
    f_nodes = f(nodes)          # f(k/n)
    BnX = bernsteinficx(f_nodes, n, X)
    approxs[n] = BnX
    sup_err = np.max(np.abs(BnX - fX))
    sup_errors.append(sup_err)
    print(f"n = {n:3d} sup-norm error = {sup_err:.6e}")
```

```python
# ---------- Figure 1: target function and selected approximations ----------
plt.figure(figsize=(10, 6))
plt plot(X, fX, label="Target function f(x)", linewidth=2, color='black')

for n in degrees:
    plt plot(X, approxs[n], label=f"Bernstein B_{n}(f)", linewidth=1)
    plt title("Stone Weierstrass (Weierstrass) demo: Bernstein polynomial approximations")
    plt xlabel("x in [0,1]")
    plt ylabel("function value")
    plt legend()
    plt grid(alpha=0.4)
    plt tight layout()
    plt show()
```

```python
# ---------- Figure 2: zoom near the kink to show uniform approximation ----------
zoom_center = 0.3
zoomRadius = 0.12
mask = (X >= zoom_center - zoom_range) & (X <= zoom_center + zoom_range)

fig = plt figure(figsize=(10, 4.5))
fig plot(X[mask], fX[mask], label="Target f(x)", linewidth=2, color='black')

for n in degrees:
    fig plot(X[mask], approxs[n][mask], label=f"B_{n}", linewidth=1)
    fig title("Zoom near the kink convergence is uniform")
    fig xlabel("x")
    fig ylabel("value")
    fig legend()
```
```python
81 plt.grid(alpha=0.4)
82 plt.tight_layout()
83 plt.show()
84
85 # ---------- Figure 3: sup-norm error vs n (log-scale) ----------
86 plt.figure(figsize=(7,5))
87 plt.plot(degrees, sup_errors, marker='o', linestyle='-')
88 plt.xscale('log')
89 plt.yscale('log')
90 plt.title("Sup-norm error (log-log scale)")
91 plt.xlabel("Degree n (log scale)")
92 plt.ylabel("Sup-norm error (log scale)")
93 plt.grid(which='both', alpha=0.4)
94 plt.tight_layout()
95 plt.show()
```

## 2.3 Depth vs. Width: Capacity Analysis

### 2.3.1 Bounding the Expressive Power

The Kolmogorov-Arnold Superposition Theorem is a foundational result in the mathematical analysis of multivariate continuous functions and their decompositions, providing a framework that underpins the expressive power of neural networks. It asserts that any continuous multivariate function can be expressed as a finite composition of continuous univariate functions and addition. It was first conjectured by Andrey Kolmogorov in 1956 and later rigorously proved by Vladimir Arnold in 1957. Formally, the theorem guarantees that any continuous multivariate function $f: [0, 1]^n \to \mathbb{R}$ can be represented as a finite composition of continuous univariate functions $\Phi_q$ and $\psi_{pq}$. Specifically, for $f(x_1, x_2, \dots, x_n)$, there exist functions $\Phi_q: \mathbb{R} \to \mathbb{R}$ and $\psi_{pq}: \mathbb{R} \to \mathbb{R}$, such that

$$f(x_1, x_2, \dots, x_n) = \sum_{q=0}^{2n} \Phi_q \left( \sum_{p=1}^{n} \psi_{pq}(x_p) \right), \qquad (2.21)$$

where the functions $\psi_{pq}(x_p)$ encode the univariate projections of the input variables $x_p$, and the outer functions $\Phi_q$ aggregate these projections into the final output. This decomposition highlights a fundamental property of multivariate continuous functions: their expressiveness can be captured through hierarchical compositions of simpler, univariate components.

#### 2.3.1.1 Literature Review of Kolmogorov-Arnold Superposition Theorem

There are some Classical References on the Kolmogorov-Arnold Superposition Theorem (KST). Kolmogorov (1957) [84] in his Foundational Paper on KST established that any continuous function of several variables can be represented as a superposition of continuous functions of a single variable and addition. This was groundbreaking because it provided a universal function decomposition method, independent of inner-product spaces. He proved that there exist functions $\phi_q$ and $\psi_q$ such that any function $f(x_1, x_2, \dots, x_n)$ can be expressed as:

$$f(x_1, \dots, x_n) = \sum_{q=1}^{2n+1} \phi_q \left( \sum_{p=1}^{n} \psi_{qp}(x_p) \right) \qquad (2.22)$$

where the $\psi_{qp}$ are univariate functions. Kolmogorov provided a mathematical basis for approximation theory and neural networks, influencing modern machine learning architectures. Arnold (1963) [86] refined Kolmogorov's theorem by proving that one can restrict the superposition to functions
of at most two variables instead of one. Arnold's formulation led to the **Kolmogorov-Arnold representation**:

$$f(x_1, \dots, x_n) = \sum_{q=1}^{2n+1} \phi_q \left( x_q + \sum_{p=1}^{n} \psi_{qp}(x_p) \right) \qquad (2.23)$$

making the theorem more suitable for practical computations. Arnold strengthened the expressivity of neural networks, inspiring alternative function representations in high-dimensional settings. Lorentz (2008) [87] in his book discusses the significance of KST in approximation theory and constructive mathematics. He provided error estimates for approximating multivariate functions using Kolmogorov-type decompositions. He showed how KST fits within Bernstein approximation theory. He helped frame KST in the context of function approximation, bridging it to computational applications. Building on this theoretical foundation, Hornik et. al. (1989) [58] demonstrated that multilayer feedforward networks are universal approximators, meaning that neural networks with a single hidden layer can approximate any continuous function. This work bridged the gap between the Kolmogorov-Arnold theorem and practical neural network design, providing a rigorous justification for the use of deep architectures. Pinkus (1999) [62] analyzed the role of KST in multilayer perceptrons (MLPs), showing how it influences function expressibility in neural networks. He demonstrated that feedforward neural networks can approximate arbitrary functions using Kolmogorov superposition. He also provided bounds on network depth and width required for universal approximation. He played a crucial role in understanding the theoretical power of deep learning. In more recent years, Montúfar, Pascanu, Cho, and Bengio (2014) [458] explored the expressive power of deep neural networks by analyzing the number of linear regions they can represent. Their work provided a modern perspective on the Kolmogorov-Arnold theorem, showing how depth enhances the ability of networks to model complex functions. Schmidt-Hieber (2020) [459] rigorously analyzed the approximation properties of deep ReLU networks, demonstrating their efficiency in approximating high-dimensional functions and further connecting the Kolmogorov-Arnold theorem to modern deep learning practices. Yarotsky (2017) [460] complemented this by providing explicit error bounds for approximating functions using deep ReLU networks, offering insights into how depth and activation functions influence approximation accuracy. Telgarsky (2016) [461] contributed to this body of work by rigorously proving that deeper networks can represent functions more efficiently than shallow ones, aligning with the hierarchical decomposition suggested by the Kolmogorov-Arnold theorem. This work provided theoretical insights into why depth is crucial in modern neural networks. Ghosh and Bhattacharya (2022) [328] propose a nested hierarchy of second-order upper bounds, $B_m$, on the system failure probability of coherent binary systems, showing that the classical Kounias-Vanmarcke-Hunter-Ditlevsen (KVHD) bound is the weakest member of this hierarchy ($m = 1$); they rigorously prove the monotonic non-increasing nature of $B_m$ with respect to $m$ and establish conditions for its strict improvement, providing a more accurate and conservative framework for reliability-based decision-making. Lu et. al. (2017) [462] explored the expressive power of neural networks from the perspective of width rather than depth, showing how width can also play a critical role in function approximation. This complemented the Kolmogorov-Arnold theorem by offering a more nuanced understanding of network design. Finally, Zhang et. al. (2021) [463] provided a rigorous analysis of how deep learning models generalize, which is closely related to their ability to approximate complex functions. While not directly about the Kolmogorov-Arnold theorem, their work contextualized these theoretical insights within the broader framework of generalization in deep learning, offering practical implications for the design and training of neural networks.

There are several very recent contributions in the Kolmogorov-Arnold Superposition Theorem (KST) (2024–2025). Guilhoto and Perdikaris (2024) [88] explored how KST can be reformulated using deep learning architectures. They proposed Kolmogorov-Arnold Networks (KANs), a new type of neural network inspired by KST. They showed that KANs outperform traditional feedforward networks in function approximation tasks. They also provided empirical evidence of KAN efficiency
in real-world datasets. They also introduced a new paradigm in machine learning, making function decomposition more interpretable. Alhafiz, M. R. et al. (2025) [89] applied KST-based networks to turbulence modeling in fluid mechanics. They demonstrated how KANs improve predictive accuracy for Navier-Stokes turbulence models. They showed a reduction in computational complexity compared to classical turbulence models. They also developed a data-driven turbulence modeling framework leveraging KST. They advanced machine learning applications in computational fluid dynamics (CFD). Lorencin, I. et al. (2024) [90] used KST-inspired neural networks for predicting propulsion system parameters in ships. They implemented KANs to model hybrid ship propulsion (Combined Diesel-Electric and Gas - CODLAG) and demonstrated a highly accurate prediction model for propulsion efficiency. They also provided a new benchmark dataset for ship propulsion research. They extended KST applications to naval engineering & autonomous systems.

<table><thead><tr><td>Paper</td><td>Main Contribution</td><td>Impact</td></tr></thead><tbody><tr><td>Kolmogorov (1957)</td><td>Original KST theorem</td><td>Laid foundation for function decomposition</td></tr><tr><td>Arnold (1963)</td><td>Refinement using 2-variable functions</td><td>Made KST more practical for computation</td></tr><tr><td>Lorentz (2008)</td><td>KST in approximation theory</td><td>Linked KST to function approximation error</td></tr><tr><td>Pinkus (1999)</td><td>KST in neural networks</td><td>Theoretical basis for deep learning</td></tr><tr><td>Perdikaris (2024)</td><td>Deep learning reinterpretation</td><td>Proposed Kolmogorov-Arnold Networks</td></tr><tr><td>Alhafiz (2025)</td><td>KST-based turbulence modeling</td><td>Improved CFD simulations</td></tr><tr><td>Lorencin (2024)</td><td>KST in naval propulsion</td><td>Optimized ship energy efficiency</td></tr></tbody></table>

### 2.3.1.2 Analysis of Kolmogorov-Arnold Superposition Theorem

In the context of neural networks, this result establishes the theoretical universality of function approximation. A neural network with a single hidden layer approximates a function $f(x_1, x_2, \dots, x_n)$ by representing it as

$$f(x_1, x_2, \dots, x_n) \approx \sum_{i=1}^{W} a_i \sigma \left( \sum_{j=1}^{n} w_{ij} x_j + b_i \right), \quad (2.24)$$

where $W$ is the width of the hidden layer, $\sigma$ is a nonlinear activation function, $w_{ij}$ are weights, $b_i$ are biases, and $a_i$ are output weights. The expressive power of such shallow networks depends critically on the width $W$, as the universal approximation theorem ensures that $W \to \infty$ suffices to approximate any continuous function arbitrarily well. However, for a fixed approximation error $\epsilon > 0$, the required width grows exponentially with the input dimension $n$, satisfying a lower bound of

$$W \geq C \cdot \epsilon^{-n}, \quad (2.25)$$

where $C$ depends on the function's Lipschitz constant. This exponential dependence, sometimes called the "curse of dimensionality," underscores the inefficiency of shallow architectures in capturing high-dimensional dependencies. The advantage of depth becomes apparent when we consider deep neural networks, which utilize hierarchical representations. A deep network with $D$ layers and width $W$ per layer constructs a function as a composition of layer-wise transformations:

$$h^{(k)} = \sigma(W^{(k)}h^{(k-1)} + b^{(k)}), \quad h^{(0)} = x, \quad (2.26)$$

where $h^{(k)}$ denotes the output of the $k$-th layer, $W^{(k)}$ is the weight matrix, $b^{(k)}$ is the bias vector, and $\sigma$ is the nonlinear activation. The final output of the network is then given by

$$f(x) \approx h^{(D)} = \sigma(W^{(D)}h^{(D-1)} + b^{(D)}). \quad (2.27)$$

The depth $D$ of the network allows it to approximate hierarchical compositions of functions. For example, if a target function $f(x)$ has a compositional structure

$$f(x) = g_1 \circ g_2 \circ \dots \circ g_D(x), \quad (2.28)$$
<figure><img src="image_21.png" /><figcaption>Figure 2.6: Kolmogorov-Arnold Superposition Theorem (Toy Example)</figcaption></figure>

where each $g_i$ is a simple function, the depth $D$ directly corresponds to the number of nested transformations. This compositional hierarchy enables deep networks to approximate functions efficiently, achieving a reduction in the required parameter count. The approximation error $\epsilon$ for a deep network decreases polynomially with $D$, satisfying

$$ \epsilon \le O\left(\frac{1}{D^2}\right), \qquad (2.29) $$

which is exponentially more efficient than the error scaling for shallow networks. In light of the Kolmogorov-Arnold theorem, the decomposition

$$ f(x_1, x_2, \dots, x_n) = \sum_{q=0}^{2n} \Phi_q \left( \sum_{p=1}^{n} \psi_{pq}(x_p) \right) \qquad (2.30) $$

demonstrates how deep networks align naturally with the structure of multivariate functions. The inner functions $\psi_{pq}$ capture local dependencies, while the outer functions $\Phi_q$ aggregate these into a global representation. This layered decomposition mirrors the depth-based structure of neural networks, where each layer learns a specific aspect of the function’s complexity. Finally, the parameter count in a deep network with $D$ layers and width $W$ per layer is given by

$$ P \le O(D \cdot W^2), \qquad (2.31) $$

whereas a shallow network requires

$$ P \ge O(W^n) \qquad (2.32) $$

parameters for the same approximation accuracy. This exponential difference in parameter count illustrates the superior efficiency of deep architectures, particularly for high-dimensional functions. By leveraging the hierarchical decomposition inherent in the Kolmogorov-Arnold theorem, deep networks achieve expressive power that scales favorably with both dimension and complexity.

### 2.3.1.3 Python Code to generate Figure 2.6

The Python code below produces the Figure 2.6 illustrating the Kolmogorov-Arnold Superposition Theorem.

```python
# """ Illustration of the Kolmogorov-Arnold Superposition Theorem: Showing how a 2D continuous function can be approximated using sums of univariate functions of linear combinations.
# """

import numpy as np
```
```python
import matplotlib.pyplot as plt

# --- Target function f(x,y) ---
def f(x, y):
    return np.sin(np.pi * x) * np.cos(np.pi * y)

# --- Toy superposition decomposition ---
# We'll handcraft psi and Phi to approximate f(x,y).
def psi11(x): return x
 def psi21(y): return y
 def psi12(x): return -x
 def psi22(y): return y

def Phi1(z): return np.sin(np.pi * z)  # expects x+y structure
 def Phi2(z): return np.cos(np.pi * z)  # expects -x+y structure

def fsuperposed(x, y):
    return 0.5*Phi1(psi11(x) + psi21(y)) + 0.5*Phi2(psi12(x) + psi22(y))

# --- Grid for plotting ---
N = 200
X = np.linspace(0, 1, N)
Y = np.linspace(0, 1, N)
XX, YY = np.meshgrid(X, Y)

F = f(XX, YY)
F_sup = f_supernposed(XX, YY)
Error = np.abs(F - F_sup)

# --- Plot results ---
fig, axes = plt.subplots(1, 3, figsize=(15, 4))

axs[0].contourf(XX, YY, F, levels=50, cmap="viridis")
axs[0].set_title("Original f(x,y)")
fig.colorbar(axs[0])

axs[1].contourf(XX, YY, F_sup, levels=50, cmap="viridis")
axs[1].set_title("Superposed Approximation")
fig.colorbar(axs[1])

axs[2].contourf(XX, YY, Error, levels=50, cmap="inferno")
axs[2].set_title("Absolute Error")
fig.colorbar(axs[2])

for ax in axes:
    ax.set_xlabel("x")
    ax.set_ylabel("y")

fig.suptitle("Kolmogorov-Arnold Superposition Theorem (Toy Example)",
             fontsize=14)

fig.tight_layout()

fig.show()
```
## 2.3.2 Fourier Analysis of Expressivity

### 2.3.2.1 Literature Review of Fourier Analysis of Expressivity

Juárez-Osorio et. al. (2024) [227] applied Fourier analysis to design quantum convolutional neural networks (QCNNs) for time series forecasting. The Fourier series decomposition helps analyze and optimize expressivity in quantum architectures, making QCNNs better at capturing periodic and non-periodic structures in data. Umeano and Kyriienko (2024) [228] introduced Fourier-based quantum feature maps that transform classical data into quantum states with enhanced expressivity. The Fourier transform plays a central role in mapping high-dimensional data efficiently while maintaining interpretability. Liu et. al. (2024) [229] extended Graph Convolutional Networks (GCNs) by integrating Fourier analysis and spectral wavelets to improve graph expressivity. It bridges the gap between frequency-domain analysis and graph embeddings, making GCNs more effective for complex data structures. Vlasic (2024) [230] presented a Fourier series-inspired feature mapping technique to encode classical data into quantum circuits. It demonstrates how Fourier coefficients can enhance the representational capacity of quantum models, leading to better compression and generalization. Kim et. al. (2024) [231] introduced Neural Fourier Modelling (NFM), a novel approach to representing time-series data compactly while preserving its expressivity. It outperforms traditional models like Short-Time Fourier Transform (STFT) in retaining long-term dependencies. Xie et. al. (2024) [232] explored how Fourier basis functions can be used to enhance the expressivity of tensor networks while maintaining computational efficiency. It establishes trade-offs between expressivity and model complexity in machine learning architectures. Liu et. al. (2024) [233] integrated spectral modulation and Fourier transforms into implicit neural representations for text-to-image synthesis. Fourier analysis improves global coherence while preserving local expressivity in generative models. Zhang (2024) [234] demonstrated how Fourier and Lock-in spectrum techniques can represent long-term variations in mechanical signals. The Fourier-based decomposition allows for more expressive representations of mechanical failures and degradation. Hamed and Lachiri (2024) [235] applied Fourier transformations to speech synthesis models, improving their ability to transfer expressive content from text to speech. Fourier series allows capturing prosody, rhythm, and tone variations effectively. Ghosh (2024) [989] presents a mathematically rigorous framework combining microstructural mechanics and continuum theory to analyze creep deformation, and derives constitutive models that integrate diffusion kinetics, dislocation dynamics, and grain-boundary behavior to predict time-dependent material deformation across diverse conditions. Ghosh (2023) [954] establishes the mathematical equivalence between the Galerkin Finite Element Method (FEM) formulation and the Central Difference formulation for the advection-diffusion equation, thereby providing a unified theoretical framework connecting finite element and finite difference time integration schemes. Ghosh (2023) [258] performs a detailed stability analysis of the second- and fourth-order Runge-Kutta methods, deriving and visualizing the corresponding stability regions in the complex $\lambda\Delta t$ plane, thereby providing deeper insight into the numerical stability characteristics of these widely used time integration schemes. Jahangiri (2022) [805] generalizes the classical open-top box optimization problem from rectangular bases to arbitrary polygons by identifying optimal corner cuts for regular and irregular prisms, and develops a new vertex-identification argument not found in prior literature. Lehmann et. al. (2024) [236] integrated Fourier-based deep learning models for seismic activity prediction. It explores the expressivity of Fourier Neural Operators (FNOs) in capturing wave propagations in different geological environments.

### 2.3.2.2 Analysis of Fourier Analysis of Expressivity

The Fourier analysis of expressivity in neural networks seeks to rigorously quantify how neural architectures, characterized by their depth and width, can approximate functions through the decomposition of those functions into their Fourier spectra. Consider a square-integrable function
$f: \mathbb{R}^d \to \mathbb{R}$, for which the Fourier transform is defined as

$$ \hat{f}(\xi) = \int_{\mathbb{R}^d} f(\mathbf{x})e^{-i2\pi\xi \cdot \mathbf{x}} d\mathbf{x} \qquad (2.33) $$

where $\xi \in \mathbb{R}^d$ represents the frequency. The inverse Fourier transform reconstructs the function as

$$ f(\mathbf{x}) = \int_{\mathbb{R}^d} \hat{f}(\xi)e^{i2\pi\xi \cdot \mathbf{x}} d\xi \qquad (2.34) $$

The magnitude $|\hat{f}(\xi)|$ reflects the energy contribution of the frequency $\xi$ to $f$. Neural networks

<figure><img src="image_22.png" /><figcaption>Figure 2.7: Target vs NN approximation (real space)</figcaption></figure>

approximate $f$ by capturing its Fourier spectrum, but the architecture fundamentally governs how efficiently this approximation can be achieved, especially in the presence of high-frequency components. For shallow networks with one hidden layer and a finite number of neurons, the universal approximation theorem establishes that

$$ f(\mathbf{x}) \approx \sum_{i=1}^{n} a_i \phi(\mathbf{w}_i \cdot \mathbf{x} + b_i) \qquad (2.35) $$

where $\phi$ is the activation function, $\mathbf{w}_i \in \mathbb{R}^d$ are weights, $b_i \in \mathbb{R}$ are biases, and $a_i \in \mathbb{R}$ are coefficients. The Fourier transform of this representation can be expressed as

$$ \hat{f}(\xi) \approx \sum_{i=1}^{n} a_i \hat{\phi}(\xi) e^{-i2\pi\xi \cdot \mathbf{b}_i} \qquad (2.36) $$

where $\hat{\phi}(\xi)$ denotes the Fourier transform of the activation function. For smooth activation functions like sigmoid or tanh, $\hat{\phi}(\xi)$ decays exponentially as $||\xi|| \to \infty$, limiting the network's ability
<figure><img src="image_23.png" /><figcaption>Figure 2.8: Fourier Spectrum: Expressivity Analysis</figcaption></figure>

to approximate functions with high-frequency content unless the width $n$ is exceedingly large. Specifically, the Fourier coefficients decay as

$$|\hat{f}(\xi)| \sim e^{-\beta \|\xi\|} \qquad (2.37)$$

where $\beta > 0$ depends on the smoothness of $\phi$. This restriction implies that shallow networks are biased toward low-frequency functions unless their width scales exponentially with the input dimension $d$. Deep networks, on the other hand, leverage their hierarchical structure to overcome these limitations. A deep network with $L$ layers recursively composes functions, producing an output of the form

$$f(\mathbf{x}) = \phi_L(\mathbf{W}^{(L)} \phi_{L-1}(\mathbf{W}^{(L-1)} \cdots \phi_1(\mathbf{W}^{(1)}\mathbf{x} + \mathbf{b}^{(1)} \cdots) + \mathbf{b}^{(L)}) \qquad (2.38)$$

where $\phi_l$ is the activation function at layer $l$, $\mathbf{W}^{(l)}$ are weight matrices, and $\mathbf{b}^{(l)}$ are bias vectors. The Fourier transform of this composition can be analyzed iteratively. If $\mathbf{h}^{(l)} = \phi_l(\mathbf{W}^{(l)}\mathbf{h}^{(l-1)} + \mathbf{b}^{(l)})$ represents the output of the $l$-th layer, then

$$\widehat{\mathbf{h}}^{(l)}(\xi) = \hat{\phi}_l(\xi) * \widehat{\mathbf{W}}^{(l)} \widehat{\mathbf{h}}^{(l-1)}(\xi) \qquad (2.39)$$

where $*$ denotes convolution and $\hat{\phi}_l$ is the Fourier transform of the activation function. The recursive application of this convolution amplifies high-frequency components, enabling deep networks to approximate functions whose Fourier spectra exhibit polynomial decay. Specifically, the Fourier coefficients of a deep network decay as

$$|\hat{f}(\xi)| \sim \|\xi\|^{-\alpha L} \qquad (2.40)$$

where $\alpha$ depends on the activation function. This is in stark contrast to the exponential decay observed in shallow networks.
The activation function plays a pivotal role in shaping the Fourier spectrum of neural networks. For example, the rectified linear unit (ReLU) $\phi(x) = \max(0, x)$ introduces significant high-frequency components into the network. The Fourier transform of the ReLU activation is given by

$$ \hat{\phi}(\xi) = \frac{1}{2\pi i \xi} \qquad (2.41) $$

which decays more slowly than the Fourier transforms of smooth activations. Consequently, ReLU-based networks are particularly effective at approximating functions with oscillatory behavior. To illustrate, consider the function

$$ f(\mathbf{x}) = \sin(2\pi\xi \cdot \mathbf{x}) \qquad (2.42) $$

A shallow network requires an exponentially large number of neurons to approximate $f$ when $\|\xi\|$ is large, but a deep network can achieve the same approximation with polynomially fewer parameters by leveraging its hierarchical structure. The expressivity of deep networks can be further quantified

<figure><img src="image_24.png" /><figcaption>Figure 2.9: Target vs Approximations (Real Space)</figcaption></figure>

by considering their ability to approximate bandlimited functions, i.e., functions $f$ whose Fourier spectra are supported on $\|\xi\| \le \omega_{max}$. For a shallow network with width $n$, the required number of neurons scales as

$$ n \sim (\omega_{\max})^d \qquad (2.43) $$

where $d$ is the input dimension. In contrast, for a deep network with depth $L$, the width scales as

$$ n \sim (\omega_{\max})^{d/L} \qquad (2.44) $$

reflecting the exponential efficiency of depth in distributing the approximation of frequency components across layers. For example, if $f(\mathbf{x}) = \cos(2\pi\xi \cdot \mathbf{x})$ with $\|\xi\| = \omega_{max}$, a deep network requires significantly fewer parameters than a shallow network to approximate $f$ to the same accuracy. A neural network with an input $\mathbf{x} \in \mathbb{R}^d$ and output $f(\mathbf{x})$ can be expressed as a sum of nonlinearly transformed weighted combinations of the input. Mathematically, this can be written as

$$ f(\mathbf{x}) = \sum_{i=1}^{m} a_i \sigma(\mathbf{w}_i^T \mathbf{x} + b_i), \qquad (2.45) $$
<figure><img src="image_25.png" /><figcaption>Figure 2.10: Fourier Spectrum: Expressivity of Architectures</figcaption></figure>

where $\sigma(\cdot)$ is the activation function, $\mathbf{w}_i \in \mathbb{R}^d$ are the weight vectors, $b_i \in \mathbb{R}$ are the biases, and $a_i$ are the output weights. To understand the spectral properties of $f(\mathbf{x})$, it is necessary to analyze the Fourier transform of the activation function $\sigma(\mathbf{w}_i^T \mathbf{x} + b_i)$, as the network output consists of a superposition of such transformed functions. The Fourier transform of $f(\mathbf{x})$ is given by

$$ \hat{f}(\mathbf{k}) = \int_{\mathbb{R}^d} f(\mathbf{x}) e^{-2\pi i \mathbf{k} \cdot \mathbf{x}} d\mathbf{x}. \qquad (2.46) $$

Substituting the definition of $f(\mathbf{x})$, we obtain

$$ \hat{f}(\mathbf{k}) = \sum_{i=1}^{m} a_i \int_{\mathbb{R}^d} \sigma(\mathbf{w}_i^T \mathbf{x} + b_i) e^{-2\pi i \mathbf{k} \cdot \mathbf{x}} d\mathbf{x}. \qquad (2.47) $$

The key term in this expression is the Fourier transform of the activation function $\sigma(\mathbf{w}^T \mathbf{x} + b)$, which we denote as

$$ \hat{\sigma}(\mathbf{k}) = \int_{\mathbb{R}^d} \sigma(\mathbf{w}^T \mathbf{x} + b) e^{-2\pi i \mathbf{k} \cdot \mathbf{x}} d\mathbf{x}. \qquad (2.48) $$

The decay rate of $\hat{\sigma}(\mathbf{k})$ determines the extent to which different frequency components are retained in the Fourier spectrum of the neural network. If $\sigma(x)$ is a smooth function, then $\hat{\sigma}(\mathbf{k})$ decays rapidly for large $\|\mathbf{k}\|$, implying that high-frequency components are suppressed. If $\sigma(x)$ is piecewise continuous or non-differentiable at some points, the decay rate of $\hat{\sigma}(\mathbf{k})$ is slower, allowing the network to capture higher frequencies. To quantify the decay properties of the Fourier transform of an activation function, consider the case of the ReLU activation function, defined as

$$ \sigma(x) = \max(0, x). \qquad (2.49) $$

Its Fourier transform is given by

$$ \hat{\sigma}(k) = \frac{1}{(2\pi i k)^2}, \qquad (2.50) $$

which exhibits a power-law decay. In contrast, the sigmoid activation function

$$ \sigma(x) = \frac{1}{1 + e^{-x}} \qquad (2.51) $$
has a Fourier transform that decays exponentially,

$$ \hat{\sigma}(k) \sim e^{-\pi|k|}. \qquad (2.52) $$

This strong decay implies that networks with sigmoid activation functions are biased toward learning only low-frequency components. Similarly, the hyperbolic tangent activation function tanh($x$) also exhibits an exponential spectral decay. For a more general class of activation functions, the decay rate of $\hat{\sigma}(k)$ can be estimated as

$$ |\hat{\sigma}(k)| \le C e^{-\alpha|k|} \qquad (2.53) $$

for some constants $C, \alpha > 0$, implying that the function is a strong low-pass filter. If $\sigma(x)$ is piece-wise smooth but non-differentiable at certain points (such as ReLU), then the Fourier transform satisfies a power-law decay

$$ |\hat{\sigma}(k)| \le C |k|^{-p} \qquad (2.54) $$

for some $p > 0$. A particularly interesting case arises when the activation function is sinusoidal, such as $\sigma(x) = \sin(x)$, for which the Fourier transform does not decay at all. The implications of these spectral properties become evident in the training dynamics of neural networks. When using gradient descent, the evolution of the Fourier coefficients of $f(\mathbf{x})$ over time follows

$$ \frac{d\hat{f}(\mathbf{k}, t)}{dt} = -\lambda(\mathbf{k})\hat{f}(\mathbf{k}, t), \qquad (2.55) $$

where $\lambda(\mathbf{k})$ is an effective learning rate for each frequency. The decay behavior of $\hat{\sigma}(\mathbf{k})$ influences $\lambda(\mathbf{k})$, meaning that activation functions with strong spectral decay impose a bottleneck on the learning of high-frequency components. From a function approximation perspective, the spectral characteristics of the activation function determine the types of functions that can be efficiently represented by a neural network. Smooth activation functions lead to approximations that are predominantly low-frequency, whereas activation functions with slower Fourier decay allow the network to approximate functions with higher-frequency content. Thus, the activation function fundamentally shapes the Fourier spectrum of neural networks, controlling the network's ability to represent and learn different frequency components.

<table><thead><tr><th>Activation Function</th><th>Fourier Decay Rate</th><th>Effect on Frequency Learning</th></tr></thead><tbody><tr><td>Sigmoid σ(x) = 1/(1+e<sup>-x</sup>)</td><td>e<sup>-π|k|</sup> (Exponential)</td><td>Strong low-pass filter, retains only low frequencies</td></tr><tr><td>Tanh σ(x) = tanh(x)</td><td>e<sup>-π|k|</sup> (Exponential)</td><td>Strong low-pass filter, smooth approximates high frequencies</td></tr><tr><td>ReLU σ(x) = max(0, x)</td><td>|k|<sup>-2</sup> (Power-law)</td><td>Allows moderate frequency learning</td></tr><tr><td>Leaky ReLU σ(x) = max(αx, x)</td><td>|k|<sup>-2</sup> (Power-law)</td><td>Similar to ReLU with slightly improved high-frequency learning</td></tr><tr><td>Sinusoidal σ(x) = sin(x)</td><td>No decay</td><td>Captures all frequencies, highly oscillatory</td></tr></tbody></table>

This table summarizes the spectral characteristics of various activation functions and their impact on frequency learning in neural networks.

In summary, the Fourier analysis of expressivity rigorously demonstrates the superiority of deep networks over shallow ones in approximating complex functions. Depth introduces a hierarchical compositional structure that enables the efficient representation of high-frequency components, while width provides a rich basis for approximating the function's Fourier spectrum. Together, these properties explain the remarkable capacity of deep neural networks to approximate functions with intricate spectral structures, offering a mathematically rigorous foundation for understanding their expressivity.
### 2.3.2.3 Python Code to generate Figures 2.7 and 2.8

The Python code below produces the Figures 2.7 and 2.8 illustrating Target vs NN approximation (real space) and the Fourier Spectrum: Expressivity Analysis.

```python
#"
Fourier analysis of expressivity: Compare Fourier spectrum of a target
function vs a shallow neural network approximation.
#"
```

```python
import numpy as np
import matplotlib.pyplot as plt
from sklearn.neural_network import MLPRegressor

# --- Target function ---
def f(x):
    # mixture of low + high frequency components
    return np.sin(2*np.pi*3*x) + 0.5*np.sin(2*np.pi*15*x)

# --- Sample training data ---
N = 512
X = np.linspace(0, 1, N).reshape(-1, 1)
y = f(X).ravel()
```

```python
# --- Train shallow neural network ---
nn = MLPRegressor(hiddenLayerSizes=(20,), activation='relu',
                    solver='lbfgs', max_iter=2000, random_state=0)
nn.fit(X, y)
y_pred = nn Predict(X)
```

```python
# --- Fourier analysis ---
def fourier_spectrumonsignal):
    # Discrete Fourier Transform
    N = len(signal)
    fft vals = np .fft .fft (signal) / N
    freqs = np .fft .fftfreq(N, d =1 /N) # normalized to [0,1]
    return freqs [:N //2], np .abs(fft vals [:N //2])

freqs, F_true = fourier_spectrum(y)
_, F_pred = fourier_spectrum(y_pred)
```

```python
# --- Plots ---
# plt .figure (figsize=(12,5))
```

```python
# Real-space comparison
plt .figure (figsize=(7,5));
```

```python
# --- Fourier spectrum comparison ---
```

```python
```
```python
51 plt.figure(figsize=(7,5));
52 plt stems(freqs, F_true, linefmt='b-', markerfmt='bo', basefmt=" ", label="Target spectrum")
53 plt stems(freqs, F_pred, linefmt='r--', markerfmt='ro', basefmt=" ", label="NN spectrum")
54 plt title("Fourier Spectrum: Expressivity Analysis")
55 plt ylabel("frequency index k")
56 plt stem(freqs, F_true, linefmt='r--', markerfmt='ro', basefmt=" ", label="NN spectrum")
57 plt legend()
58 plt grid(alpha=0.4)
59 plt show()
```

### 2.3.2.4 Python Code to generate Figures 2.9 and 2.10

The Python code below produces the Figures 2.9 and 2.10 illustrating the Target vs Approximations (Real Space) and Fourier Spectrum: Expressivity of Architectures respectively.

```python
1 """
2 Compare Fourier expressivity of different neural network architectures ( shallow vs deep, ReLU vs tanh).
3 """
4
5 import numpy as np
6 import matplotlib.pyplot as plt
7 fromsklearn.neural_networkimportMLPRegressor
8
9 # --- Target function with multiple frequencies ---
10 def f(x):
11     return np.sin(2* np.pi*3*x) + 0.5* np.sin(2* np.pi*15*x)
12
13 # --- Training data ---
14 N = 512
15 X = np.linspace(0, 1, N).reshape(-1, 1)
16 Y = f(X).ravel()
17
18 # --- Helper: compute Fourier spectrum ---
19 def fourier_spectrum是有
20     N = len是有
21     fft vals = np.fft.fft是有
22     freqs = np.fft.fftfreq是有
23     return freqs[:N//2], np. abs(fft vals[:N//2])
24
25 freqs, F_true = fourier_spectrum(Y)
26
27 # --- Define architectures ---
28 architectures = {
29     "Shallow ReLU (1x20)" : MLPRe regressor(hidden layer sizes=(20,),
30                                         activation='relu',
31                                         solver='lbfgs', max_iter=3000,
32                                         random_state=0),
33     "Deep ReLU (3x20)" : MLPRe regressor(hidden layer sizes=(20,20,20),
34                                         activation='relu',
35                                         solver='lbfgs', max_iter=3000,
36                                         random_state=0),
```
```python
33 "Shallow tanh (1x20)": MLPRegressor(hiddenLayerSizes=(20, ),
34                                     activation='tanh',
35                                     solver='lbfgs', max_iter=3000,
36                                     random_state=0),
36 "Deep tanh (3x20)": MLPRegressor(hiddenLayerSizes=(20, 20, 20),
37                                     activation='tanh',
38                                     solver='lbfgs', max_iter=3000,
39                                     random_state=0),
40 # --- Train and collect spectra ---
41 spectra = {}
42 preds = {}
43 for name, model in architectures.items():
44     model.fit(X, y)
45     y_pred = model Predict(X)
46     preds[name] = y_pred
47     F_pred = fourierSpectrum(y_pred)
48     spectra[name] = F_pred
49     print(f"{name} trained.")
50
51 # --- Plot: function fits ---
52 plt.figure(figsize=(12, 6))
53 plt.Plot(X, y, 'k', lw=2, label="Target f(x)")
54 for name, y_pred in preds.items():
55     plt.Plot(X, y_pred, label=name, lw=1.5)
56     plt.TITLE("Target vs Approximations (Real Space)")
57 plt.Xlabel("x")
58 plt.Ylabel("f(x)")
59 plt.Legend()
60 plt GRID(alpha=0.4)
61 plt.TightLayout()
62 plt.Show()
63
64 # --- Plot: Fourier spectra comparison ---
65 plt.Plot(figsize=(12, 6))
66 plt.Plot(freqs, F_true, 'k-', lw=2, label="Target Spectrum")
67 for name, F_pred in spectra.items():
68     plt.Plot(freqs, F_pred, label=name)
69     plt.TITLE("Fourier Spectrum: Expressivity of Architectures")
70 plt.Xlabel("frequency index k")
71 plt.Ylabel("Amplitude")
72 plt.Legend()
73 plt GRID(alpha=0.4)
74 plt.TightLayout()
75 plt.Show()
```

### 2.3.3 Fourier transforms of various activation functions

To rigorously derive the Fourier transforms of various activation functions, let us consider the Fourier transform definition:

$$ \hat{\sigma}(k) = \int_{-\infty}^{\infty} \sigma(x)e^{-2\pi ikx} dx. \qquad (2.56) $$
We will explicitly derive $\hat{\sigma}(k)$ for each activation function in a mathematically rigorous manner. This will include the Sigmoid, Tanh, ReLU, Leaky ReLU, and Sinusoidal activation functions.

### 2.3.3.1 Fourier Transform of the Sigmoid Function

The sigmoid activation function is given by:

$$ \sigma(x) = \frac{1}{1 + e^{-x}}. \qquad (2.57) $$

The Fourier transform of $\sigma(x)$ is:

$$ \hat{\sigma}(k) = \int_{-\infty}^{\infty} \frac{1}{1+e^{-x}} e^{-2\pi ikx} dx. \qquad (2.58) $$

Rewriting $\sigma(x)$ using the identity:

<figure><img src="image_26.png" /><figcaption>Figure 2.11: Sigmoid Activation Function: $\sigma(x) = \frac{1}{1+e^{-x}}$</figcaption></figure>

$$ \sigma(x) = \frac{e^x}{1+e^x}, \qquad (2.59) $$

the Fourier integral becomes:

$$ \hat{\sigma}(k) = \int_{-\infty}^{\infty} \frac{e^x}{1+e^x} e^{-2\pi ikx} dx. \qquad (2.60) $$

Setting $u = e^x$, so that $du = e^x dx$, we transform the integral:

$$ \hat{\sigma}(k) = \int_{0}^{\infty} \frac{u}{1+u} u^{-2\pi ik} \frac{du}{u}. \qquad (2.61) $$

Simplifying:

$$ \hat{\sigma}(k) = \int_{0}^{\infty} \frac{u^{-2\pi ik+1}}{1+u} du. \qquad (2.62) $$

This integral is well-known and can be computed using contour integration, yielding:

$$ \hat{\sigma}(k) = \pi \operatorname{csch}(\pi k), \qquad (2.63) $$

where $\operatorname{csch}(x) = \frac{1}{\sinh(x)}$ is the hyperbolic cosecant function. For large $|k|$, this decays as:

$$ \hat{\sigma}(k) \sim e^{-\pi|k|}. \qquad (2.64) $$
### 2.3.3.2 Fourier Transform of the Hyperbolic Tangent Function

The tanh activation function is:

$$ \sigma(x) = \tanh(x) = \frac{e^x - e^{-x}}{e^x + e^{-x}} \qquad (2.65) $$

The Fourier transform is:

<figure><img src="image_27.png" /><figcaption>Figure 2.12: Hyperbolic Tangent Activation Function: $\sigma(x) = \frac{e^x - e^{-x}}{e^x + e^{-x}}$</figcaption></figure>

$$ \hat{\sigma}(k) = \int_{-\infty}^{\infty} \frac{e^x - e^{-x}}{e^x + e^{-x}} e^{-2\pi i kx} dx. \qquad (2.66) $$

Using the identity:

$$ \tanh(x) = 1 - 2\sigma(-x), \qquad (2.67) $$

and leveraging the result from the sigmoid case, we obtain:

$$ \hat{\sigma}(k) = -2\pi i \operatorname{csch}(\pi k). \qquad (2.68) $$

Again, for large $|k|$:

$$ \hat{\sigma}(k) \sim e^{-\pi|k|}. \qquad (2.69) $$

### 2.3.3.3 Fourier Transform of the ReLU Function

The ReLU activation function is:

$$ \sigma(x) = \max(0, x) \qquad (2.70) $$

Thus, the Fourier transform is:

$$ \hat{\sigma}(k) = \int_{-\infty}^{\infty} xH(x)e^{-2\pi i kx} dx, \qquad (2.71) $$

where $H(x)$ is the Heaviside step function. Since $H(x)$ is zero for $x < 0$, the integral simplifies to:

$$ \hat{\sigma}(k) = \int_{0}^{\infty} xe^{-2\pi i kx} dx. \qquad (2.72) $$
<figure><img src="image_28.png" /><figcaption>Figure 2.13: ReLU Activation Function: $\sigma(x) = \max(0, x)$</figcaption></figure>

Using integration by parts with $u = x$ and $dv = e^{-2\pi ikx}dx$, we find:

$$ \hat{\sigma}(k) = \frac{1}{(2\pi ik)^2}. \qquad (2.73) $$

This exhibits a power-law decay:

$$ \hat{\sigma}(k) \sim |k|^{-2}. \qquad (2.74) $$

### 2.3.3.4 Fourier Transform of the Leaky ReLU Function

The Leaky ReLU function is:

$$ \sigma(x) = \begin{cases} x, & x \ge 0, \\ \alpha x, & x < 0 \end{cases} \qquad (2.75) $$

The Fourier transform is:

$$ \hat{\sigma}(k) = \int_{-\infty}^{\infty} \sigma(x)e^{-2\pi ikx}dx. \qquad (2.76) $$

Splitting the integral:

$$ \hat{\sigma}(k) = \int_{0}^{\infty} xe^{-2\pi ikx}dx + \alpha \int_{-\infty}^{0} xe^{-2\pi ikx}dx. \qquad (2.77) $$

Computing these integrals using integration by parts gives:

$$ \hat{\sigma}(k) = \frac{1 + \alpha}{(2\pi ik)^2}. \qquad (2.78) $$

Thus, the decay is still $|k|^{-2}$, but the amplitude is modulated by $(1 + \alpha)$.

### 2.3.3.5 Fourier Transform of the Sinusoidal Activation Function

For $\sigma(x) = \sin(x)$, the Fourier transform is:

$$ \hat{\sigma}(k) = \int_{-\infty}^{\infty} \sin(x)e^{-2\pi ikx}dx. \qquad (2.79) $$
<figure><img src="image_29.png" /><figcaption>Figure 2.14: Leaky ReLU Activation Function</figcaption></figure>

Using the Fourier transform properties of sine:

$$ \hat{\sigma}(k) = \frac{i}{2} \left[ \delta(k - \frac{1}{2\pi}) - \delta(k + \frac{1}{2\pi}) \right]. \qquad (2.80) $$

This shows that the sinusoidal activation function contains only two specific frequencies and does not decay, meaning it retains high-frequency components perfectly.

These rigorous derivations confirm how different activation functions influence the Fourier spectrum of neural networks. The decay properties play a crucial role in determining the network's ability to learn high- or low-frequency functions.

## 2.4 The Connection Between Different Mathematics Problems and Deep Learning

### 2.4.1 Basel Problem and Deep Learning

The **Basel Problem** is a famous question in mathematical analysis that asks for the sum of the reciprocals of the squares of natural numbers:

$$ \sum_{n=1}^{\infty} \frac{1}{n^2}. \qquad (2.81) $$

Leonhard Euler solved this problem in the 18th century and astonishingly found that the sum converges to $\frac{\pi^2}{6}$, meaning:

$$ \sum_{n=1}^{\infty} \frac{1}{n^2} = \frac{\pi^2}{6}. \qquad (2.82) $$

This result is remarkable because it links a seemingly simple infinite sum to the fundamental mathematical constant $\pi$, which is deeply connected to trigonometry and Fourier analysis. A
<figure><img src="image_30.png" /><figcaption>Figure 2.15: Sinusoidal Activation Function $\sigma(x) = \sin(x)$</figcaption></figure>

summary of notable proofs given by mathematicians to the Basel problem is given by Ghosh (2020) [698]. A solution to the Basel problem using the Calculus of Residues is given by Ghosh (2021) [800]. In modern mathematics, the Basel Problem's solution is understood in the broader context of **zeta functions**, where the Riemann zeta function $\zeta(s)$ is defined as:

$$ \zeta(s) = \sum_{n=1}^{\infty} \frac{1}{n^s}. \qquad (2.83) $$

For $s = 2$, we recover the Basel Problem. The significance of such summations extends far beyond pure mathematics; they appear in areas such as **signal processing**, **numerical analysis**, and **machine learning**. In particular, deep learning models, which fundamentally rely on function approximation and optimization, exhibit connections to the types of series that arise in the Basel Problem, particularly in relation to Fourier series and the decay of spectral components in neural network approximations.

### 2.4.1.1 The Basel Problem, Fourier Series, and Function Approximation in Deep Learning

One of the most natural places where the Basel Problem appears in applied mathematics is in **Fourier series**. In Fourier analysis, a periodic function can be decomposed into a sum of sinusoidal functions with different frequencies. The coefficients of these sinusoidal components determine the function's structure, and their squared magnitudes are directly tied to the Basel-type sums. Specifically, for a function $f(x)$ with Fourier expansion:

$$ f(x) = \sum_{n=1}^{\infty} a_n \cos(nx) + b_n \sin(nx), \qquad (2.84) $$

Parseval's theorem states that the total energy of the function is given by:

$$ \sum_{n=1}^{\infty} (a_n^2 + b_n^2). \qquad (2.85) $$
If the function is smooth, the coefficients $a_n, b_n$ decay rapidly, often behaving like $\frac{1}{n^2}$ or faster. The Basel Problem is thus closely related to the **rate of decay of Fourier coefficients**, which directly impacts how well a function can be approximated using a truncated Fourier series. Deep learning, particularly in neural network function approximation, shares a fundamental connection with Fourier series. Neural networks with sufficient width and depth can approximate arbitrary functions, and recent research has shown that **deep networks tend to first learn lower-frequency components before higher frequencies**—a phenomenon known as **spectral bias**. The Basel-type decay of Fourier coefficients is precisely what governs the smoothness of these approximations, meaning that understanding series like the Basel Problem helps us characterize how neural networks generalize and learn complex functions.

### 2.4.1.2 The Role of the Basel Problem in Regularization and Weight Decay

Another important connection between the Basel Problem and deep learning emerges in the context of **regularization techniques**, particularly **L2 regularization**, also known as **weight decay**. Regularization techniques in neural networks help prevent overfitting by penalizing large weight magnitudes, encouraging smoother function approximations. Mathematically, L2 regularization adds a penalty term to the loss function:

$$ \mathcal{L}_{\text{reg}} = \lambda \sum_{i} w_{i}^{2}, \qquad (2.86) $$

where $w_i$ are the weights of the network, and $\lambda$ is a regularization parameter. This penalty ensures that the network does not learn excessively large weights, which could lead to high-frequency oscillations in the approximated function. The sum of squared weights in L2 regularization resembles the Basel-type summation of reciprocal squares. Just as the Basel Problem reflects the decay of Fourier coefficients for smooth functions, L2 regularization ensures that the learned function maintains smoothness by controlling the magnitude of weights. Thus, there is a deep analogy: **both involve penalizing high-frequency components to favor smooth, well-behaved solutions.**

### 2.4.1.3 Spectral Bias in Deep Learning and the Basel Problem

One fascinating empirical observation in deep learning is that neural networks naturally prioritize learning lower-frequency components of a function first. This phenomenon, known as **spectral bias**, is well-studied in the context of deep learning theory. It means that when a neural network is trained on a function, it first learns the dominant, low-frequency parts before capturing finer details at higher frequencies. This is strikingly similar to how the Basel Problem reflects the fundamental nature of **low-frequency dominance in function representations**. If we consider a function that can be expanded in a Fourier series, the Basel sum describes the contribution of lower-frequency components, which dominate the structure of the function. Similarly, in deep learning, neural networks exhibit behavior where early learning stages prioritize simple, smooth approximations, just as Fourier coefficients with lower indices contribute the most to a function's overall shape.

Thus, the decay of Fourier coefficients (as described by the Basel Problem) provides an intuition for why neural networks prefer lower frequencies. This insight is valuable in designing deep learning architectures, especially in areas like **implicit neural representations (INRs)**, where Fourier features are explicitly incorporated into the model design to control the spectral bias. While the Basel Problem originates in pure mathematics, its influence extends deeply into areas of applied mathematics, including function approximation, signal processing, and neural networks. The problem's solution reveals a fundamental truth about how series behave, and this truth manifests in Fourier series, weight decay regularization, and spectral bias in deep learning. The key takeaway is that **the Basel Problem describes the natural decay of Fourier coefficients in function**
**approximations**, and similar decay patterns emerge in deep learning when training networks to approximate complex functions. Euler's result continues to play a role in shaping our understanding of function representations in neural networks today.
## 3 Training Dynamics and NTK Linearization

Training dynamics refers to the evolution of a neural network's parameters during optimization, typically via gradient descent (GD) or stochastic gradient descent (SGD). Given a neural network $f(\mathbf{x}; \boldsymbol{\theta})$ with parameters $\boldsymbol{\theta} \in \mathbb{R}^P$, trained on a dataset $\{(\mathbf{x}_i, y_i)\}_{i=1}^N$, the loss function $\mathcal{L}(\boldsymbol{\theta})$ is minimized using the gradient flow (continuous-time approximation of GD):

$$ \frac{d\boldsymbol{\theta}(t)}{dt} = -\nabla_{\boldsymbol{\theta}}\mathcal{L}(\boldsymbol{\theta}(t)), \quad \text{where} \quad \mathcal{L}(\boldsymbol{\theta}) = \frac{1}{2}\sum_{i=1}^{N}(f(\mathbf{x}_i; \boldsymbol{\theta}) - y_i)^2. \qquad (3.1) $$

The dynamics of the network's predictions $\mathbf{f}_t = (f(\mathbf{x}_1; \boldsymbol{\theta}(t)), \dots, f(\mathbf{x}_N; \boldsymbol{\theta}(t)))$ can be described by:

$$ \frac{d\mathbf{f}_t}{dt} = -\mathbf{H}_t(\mathbf{f}_t - \mathbf{y}), \qquad (3.2) $$

where $\mathbf{H}_t$ is the Neural Tangent Kernel (NTK) matrix, defined as:

$$ \mathbf{H}_t = \nabla_{\boldsymbol{\theta}}\mathbf{f}_t \cdot (\nabla_{\boldsymbol{\theta}}\mathbf{f}_t)^T. \qquad (3.3) $$

The NTK governs how small changes in parameters affect the network's outputs, determining the trajectory of learning in function space. For wide neural networks, the NTK remains approximately

<figure><img src="image_31.png" /><figcaption>Figure 3.1: Training Dynamics vs NTK Linearization</figcaption></figure>

constant during training, leading to a simplified linearized dynamics approximation. Let $\boldsymbol{\theta}_0$ denote the initial parameters (typically randomly sampled). The first-order Taylor expansion of $f(\mathbf{x}; \boldsymbol{\theta})$ around $\boldsymbol{\theta}_0$ yields:

$$ f_{\text{lin}}(\mathbf{x}; \boldsymbol{\theta}) = f(\mathbf{x}; \boldsymbol{\theta}_0) + \nabla_{\boldsymbol{\theta}} f(\mathbf{x}; \boldsymbol{\theta}_0)^T (\boldsymbol{\theta} - \boldsymbol{\theta}_0). \qquad (3.4) $$
## 3.1. PYTHON CODE TO GENERATE FIGURE 3.1 ILLUSTRATING THE TRAINING DYNAMICS VS NTK

Under this linearization, the training dynamics simplify to kernel gradient descent with a fixed kernel $H_0$, the NTK at initialization. The evolution of predictions becomes:

$$f_t = y + e^{-H_0 t}(f_0 - y), \qquad (3.5)$$

where $f_0$ is the initial network output. This shows that the network behaves like a linear model in the infinite-width limit, with convergence governed by the eigendecomposition of $H_0$. The NTK theory provides a precise characterization of generalization and optimization in overparameterized neural networks, bridging kernel methods and deep learning. The key insight is that despite their nonlinearity, wide neural networks evolve in a regime where their training dynamics are effectively governed by a static, deterministic kernel.

## 3.1 Python Code to generate Figure 3.1 illustrating the Training Dynamics vs NTK Linearization

The Python code below produces the Figure 3.1 illustrating the Training Dynamics vs NTK Linearization.

```python
import matplotlib.pyplot as plt
import numpy as np
import torch
import torch.nn as nn
import torch.optim as optim
from torch.autograd import grad

torch.manual_seed(0)
n = 100
x = torch.linspace(-1, 1, n).reshape(-1, 1)
y_true = torch.sin(3 * x) # target function

class SimpleNet nn.Module:
    def __init__(self, hidden=50):
        super().__init__()
        self.fc1 = nn.Linear(1, hidden)
        self relu = nn.Tanh()
        self.fc2 = nn.Linear(hidden, 1)

    def forward(self, x):
        return self.fc2(self relu(self.fc1(x)))

# Initialize
net = SimpleNet()
optimiser = optim.SGD(net.parameters(), lr=0.1)

# ....................................
# 3. Train & record dynamics
# ....................................
loss_fn = nn.MSELoss()
```
epochs = 50
train_lostes = []
preds_dynamics = []

for epoch in range(epochs):
    optimizer.zeroGrad()
    y_pred = net(x)
    loss = loss_fn(y_pred, y_true)
    loss backward()
    optimizer.step()

    train_lostes.append(loss.item())
    preds_dynamics.append(y_pred detach(). numpy())

# ==============================
# 4. Compute NTK Linearization (fixed)
# ==============================

# Initial output
f0 = netisman().detach(). numpy()

# Jacobian features: J in R^{n x p}
J = jacobian_features(net man, x).detach(). numpy()

# NTK kernel
K = J @ J.T # (n x n)

# Regularization for stability
K_inv = np lINALG. Pnv(K + 1e-6 * np . eye(n))

# Alpha lives in R^n
alpha = K_inv @ (y_true . numpy() - f0)

# NTK prediction
f man = f0 + K @ alpha # matches n x 1

# ==============================
# 5. Plot Training Dynamics vs NTK
# ==============================

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

# True function
p1t . plot (x, y_true, 'k--', label="True Function")

# Final NN dynamics
p1t . plot (x, preds_dynamics [-1], 'r', label="NN (Trained)")

# NTK linearization
p1t . plot (x, f man, 'b', label="NTK Linearization")

p1t . title ("Training Dynamics vs NTK Linearization")
p1t .xlabel("x")
p1t .ylabel("y")
p1t .legend()
p1t .show()
## 3.2 Literature Review of Training Dynamics and NTK Linearization

Trevisan et. al. [92] investigated how knowledge distillation can be analyzed using the Neural Tangent Kernel (NTK) framework and demonstrated that under certain conditions, the training dynamics of a student model in knowledge distillation closely follow NTK linearization. They explored how NTK affects generalization and feature transfer in the distillation process. They provided theoretical insight into why knowledge distillation improves performance in deep networks. Bonfanti et. al. (2024) [93] studied how NTK behaves in the nonlinear regime, particularly in Physics-Informed Neural Networks (PINNs). They showed that when PINNs operate outside the NTK regime, their performance degrades due to high sensitivity to initialization and weight updates. They established conditions under which NTK linearization is insufficient for PINNs, emphasizing the need for nonlinear adaptations. They provided practical guidelines for designing PINNs that maintain stable training dynamics. Jacot et. al. (2018) [94] introduced the Neural Tangent Kernel (NTK) as a fundamental framework for analyzing infinite-width neural networks. They proved that as width approaches infinity, neural networks evolve as linear models governed by the NTK. They derived generalization bounds for infinitely wide networks and connected training dynamics to kernel methods. They established NTK as a core tool in deep learning theory, leading to further developments in training dynamics research. Lee et. al. (2019) [95] extended NTK theory to arbitrarily deep networks, showing that even deep architectures behave as linear models under gradient descent and proved that training dynamics remain stable regardless of network depth when width is sufficiently large. They explored practical implications for initializing and optimizing deep networks. They strengthened NTK theory by confirming its validity beyond shallow networks. Yang and Hu (2022) [96] challenged the conventional NTK assumption that feature learning is negligible in infinite-width networks and showed that certain activation functions can induce nontrivial feature learning even in infinite-width regimes. They suggested that feature learning can be integrated into NTK theory, opening new directions in kernel-based deep learning research. Xiang et. al. (2023) [97] investigated how finite-width effects impact training dynamics under NTK assumptions and showed that finite-width networks deviate from NTK predictions due to higher-order corrections in weight updates. They derived corrections to NTK theory for practical networks, improving its predictive power for real-world architectures. They refined NTK approximations, making them more applicable to modern deep-learning models. Lee et. al. (2019) [98] extended NTK linearization to deep convolutional networks, analyzing their training dynamics under infinite width and showed how locality and weight sharing in CNNs impact NTK behavior. They also demonstrated practical consequences for CNN training in real-world applications. They bridged NTK theory and convolutional architectures, providing new theoretical tools for CNN analysis.

## 3.3 Gradient Flow and Stationary Points

### 3.3.1 Literature Review of Gradient Flow and Stationary Points

Goodfellow et. al. (2016) [120] provided a comprehensive overview of deep learning, including a detailed discussion of gradient-based optimization methods. It rigorously explains the dynamics of gradient descent in the context of neural networks, covering topics such as backpropagation, vanishing gradients, and saddle points. The book also discusses the role of learning rates, momentum, and adaptive optimization methods in shaping the trajectory of gradient flow. Sra et. al. (2012) [492] included several chapters dedicated to the theoretical and practical aspects of gradient-based optimization in machine learning. It provides rigorous mathematical treatments of gradient flow dynamics, including convergence analysis, the impact of stochasticity in stochastic gradient de-
scent (SGD), and the geometry of loss landscapes in high-dimensional spaces. Choromanska et. al. (2015) [493] rigorously analyzed the loss surfaces of deep neural networks. It demonstrates that the loss landscape is highly non-convex but contains a large number of local minima that are close in function value to the global minimum. The paper provides insights into how gradient flow navigates these complex landscapes and why it often converges to satisfactory solutions despite the non-convexity. Arora et al. (2019) [494] provided a theoretical framework for understanding the dynamics of gradient descent in deep neural networks. It rigorously analyzes the role of over-parameterization in enabling gradient flow to converge to global minima, even in the absence of explicit regularization. The paper also explores the implicit regularization effects of gradient descent and their impact on generalization. Du et. al. (2019) [485] establishes theoretical guarantees for the convergence of gradient descent to global minima in overparameterized neural networks. It rigorously proves that gradient flow can efficiently minimize the training loss to zero, even in the presence of non-convexity, by leveraging the high-dimensional geometry of the loss landscape. The authors provided a rigorous analysis of the exponential convergence of gradient descent in overparameterized neural networks. It shows that the gradient flow dynamics are characterized by a rapid decrease in the loss function, driven by the alignment of the network's parameters with the data. The paper also discusses the role of initialization in shaping the trajectory of gradient flow. Zhang et al. (2017) [463] challenged traditional notions of generalization in deep learning. It rigorously demonstrates that deep neural networks can fit random labels, suggesting that the dynamics of gradient flow are not solely driven by the data distribution but also by the implicit biases of the optimization algorithm. The paper highlights the importance of understanding how gradient flow interacts with the architecture and initialization of neural networks. Baratin et. al. (2020) [495] explored the implicit regularization effects of gradient flow in deep learning from the perspective of function space. It rigorously demonstrates that gradient descent in overparameterized models tends to converge to solutions that minimize certain norms or complexity measures, providing insights into why these models generalize well despite their capacity to overfit. Balduzzi et al. (2018) [496] extended the analysis of gradient flow to multi-agent optimization problems, such as those encountered in generative adversarial networks (GANs). It rigorously characterizes the dynamics of gradient descent in games, highlighting the role of rotational forces and the challenges of convergence in non-cooperative settings. The paper provides tools for understanding how gradient flow behaves in complex, interactive learning scenarios. Allen-Zhu et al. (2019) [487] provided a rigorous convergence theory for deep learning models trained with gradient descent. It shows that overparameterization enables gradient flow to avoid bad local minima and converge to global minima efficiently. The paper also analyzes the role of initialization, step size, and network depth in shaping the dynamics of gradient descent.

### 3.3.2 Analysis of Gradient Flow and Stationary Points

The dynamics of gradient flow in neural network training are fundamentally governed by the continuous evolution of parameters $\theta(t)$ under the influence of the negative gradient of the loss function, expressed as

$$ \frac{d\theta(t)}{dt} = -\nabla_{\theta}\mathcal{L}(\theta(t)). \qquad (3.6) $$

The loss function, typically of the form

$$ \mathcal{L}(\theta) = \frac{1}{2n} \sum_{i=1}^{n} \|f(x_i; \theta) - y_i\|^2, \qquad (3.7) $$

measures the discrepancy between the network's predicted outputs $f(x_i; \theta)$ and the true labels $y_i$. At stationary points of the flow, the condition

$$ \nabla_{\theta}\mathcal{L}(\theta^*) = 0 \qquad (3.8) $$
holds, indicating that the gradient vanishes. To classify these stationary points, the Hessian matrix $H = \nabla_{\theta}^{2}\mathcal{L}(\theta)$ is examined. For eigenvalues $\{\lambda_{i}\}$ of $H$, the nature of the stationary point is determined: $\lambda_{i} > 0$ for all $i$ corresponds to a local minimum, $\lambda_{i} < 0$ for all $i$ to a local maximum, and mixed signs indicate a saddle point. Under gradient flow $\frac{d\theta(t)}{dt} = -\nabla_{\theta}\mathcal{L}(\theta(t))$, the trajectory converges to critical points:

$$ \lim_{t \to \infty} \| \nabla_{\theta} \mathcal{L}(\theta(t)) \| = 0. \qquad (3.9) $$

The gradient flow also governs the temporal evolution of the network's predictions $f(x; \theta(t))$. A Taylor series expansion of $f(x; \theta)$ about an initial parameter $\theta_{0}$ gives:

$$ f(x; \theta) = f(x; \theta_0) + J_f(x; \theta_0)(\theta - \theta_0) + \frac{1}{2}(\theta - \theta_0)^{\top} H_f(x; \theta_0)(\theta - \theta_0) + \mathcal{O}(\|\theta - \theta_0\|^3), \quad (3.10) $$

where $J_f(x; \theta_0) = \nabla_{\theta} f(x; \theta_0)$ is the Jacobian and $H_f(x; \theta_0)$ is the Hessian of $f(x; \theta)$ with respect to $\theta$. In the NTK (neural tangent kernel) regime, higher-order terms are negligible due to the large

<figure><img src="image_32.png" /><figcaption>Figure 3.2: Gradient Flow and Stationary Points</figcaption></figure>

parameterization of the network, and the linear approximation suffices:

$$ f(x; \theta) \approx f(x; \theta_0) + J_f(x; \theta_0)(\theta - \theta_0). \qquad (3.11) $$

Under gradient flow, the time derivative of the network's predictions is given by:

$$ \frac{df(x; \theta(t))}{dt} = J_f(x; \theta(t)) \frac{d\theta(t)}{dt}. \qquad (3.12) $$

Substituting the parameter dynamics $\frac{d\theta(t)}{dt} = -\nabla_{\theta}\mathcal{L}(\theta(t)) = -\sum_{i=1}^{n}(f(x_i; \theta(t)) - y_i)J_f(x_i; \theta(t))$, this becomes:

$$ \frac{df(x; \theta(t))}{dt} = - \sum_{i=1}^{n} J_f(x; \theta(t)) J_f(x_i; \theta(t))^{\top} (f(x_i; \theta(t)) - y_i). \qquad (3.13) $$
Defining the NTK as $\mathcal{K}(x, x'; \theta) = J_f(x; \theta)J_f(x'; \theta)^\top$, and assuming constancy of the NTK during training ($\mathcal{K}(x, x'; \theta) \approx \mathcal{K}_0(x, x')$), the evolution equation simplifies to:

$$ \frac{df(x; \theta(t))}{dt} = - \sum_{i=1}^{n} \mathcal{K}_0(x, x_i)(f(x_i; \theta(t)) - y_i). \quad (3.14) $$

Rewriting in matrix form, let $\mathbf{f}(t) = [f(x_1; \theta(t)), \dots, f(x_n; \theta(t))]^\top$ and $\mathbf{y} = [y_1, \dots, y_n]^\top$. The NTK matrix $\mathcal{K}_0 \in \mathbb{R}^{n \times n}$ evaluated at initialization defines the system:

$$ \frac{d\mathbf{f}(t)}{dt} = -\mathcal{K}_0(\mathbf{f}(t) - \mathbf{y}). \quad (3.15) $$

The solution to this linear system is:

$$ \mathbf{f}(t) = e^{-\mathcal{K}_0 t} \mathbf{f}(0) + (\mathbf{I} - e^{-\mathcal{K}_0 t})\mathbf{y}. \quad (3.16) $$

As $t \to \infty$, the predictions converge to the labels: $\mathbf{f}(t) \to \mathbf{y}$, implying zero training error. The eigenvalues of $\mathcal{K}_0$ determine the rates of convergence. Diagonalizing $\mathcal{K}_0$ as $\mathcal{K}_0 = Q\Lambda Q^\top$, where $Q$ is orthogonal and $\Lambda = \text{diag}(\lambda_1, \dots, \lambda_n)$, the dynamics in the eigenbasis are:

$$ \frac{d\tilde{\mathbf{f}}(t)}{dt} = -\Lambda(\tilde{\mathbf{f}}(t) - \tilde{\mathbf{y}}), \quad (3.17) $$

with $\tilde{\mathbf{f}}(t) = Q^\top \mathbf{f}(t)$ and $\tilde{\mathbf{y}} = Q^\top \mathbf{y}$. Solving, we obtain:

$$ \tilde{\mathbf{f}}(t) = e^{-\Lambda t} \tilde{\mathbf{f}}(0) + (\mathbf{I} - e^{-\Lambda t}) \tilde{\mathbf{y}}. \quad (3.18) $$

Each mode decays exponentially with a rate proportional to the eigenvalue $\lambda_i$. Modes with larger $\lambda_i$ converge faster, while smaller eigenvalues slow convergence.

The NTK framework thus rigorously explains the linearization of training dynamics in overparameterized neural networks. This linear behavior ensures that the optimization trajectory remains within a convex region of the parameter space, leading to both convergence and generalization. By leveraging the constancy of the NTK, the complexity of nonlinear neural networks is reduced to an analytically tractable framework that aligns closely with empirical observations.

### 3.3.3 Hessian Structure

The Hessian matrix, $H(\theta) = \nabla_{\theta}^2 \mathcal{L}(\theta)$, serves as a critical construct in the mathematical framework of optimization, capturing the second-order partial derivatives of the loss function $\mathcal{L}(\theta)$ with respect to the parameter vector $\theta \in \mathbb{R}^d$. Each element $H_{ij} = \frac{\partial^2 \mathcal{L}(\theta)}{\partial \theta_i \partial \theta_j}$ reflects the curvature of the loss surface along the $(i, j)$-direction. The symmetry of $H(\theta)$, guaranteed by the Schwarz theorem under the assumption of continuous second partial derivatives, implies $H_{ij} = H_{ji}$. This property ensures that the eigenvalues $\lambda_1, \lambda_2, \dots, \lambda_d$ of $H(\theta)$ are real and the eigenvectors $\mathbf{v}_1, \mathbf{v}_2, \dots, \mathbf{v}_d$ are orthogonal, satisfying the eigenvalue equation

$$ H(\theta)\mathbf{v}_i = \lambda_i \mathbf{v}_i \quad \text{for all } i. \quad (3.19) $$

The behavior of the loss function around a specific parameter value $\theta_0$ can be rigorously analyzed using a second-order Taylor expansion. This expansion is given by:

$$ \mathcal{L}(\theta) = \mathcal{L}(\theta_0) + (\theta - \theta_0)^{\top} \nabla_{\theta} \mathcal{L}(\theta_0) + \frac{1}{2} (\theta - \theta_0)^{\top} H(\theta_0) (\theta - \theta_0) + \mathcal{O}(\|\theta - \theta_0\|^3). \quad (3.20) $$

Here, the term $(\theta - \theta_0)^{\top} \nabla_{\theta} \mathcal{L}(\theta_0)$ represents the linear variation of the loss, while the quadratic term $\frac{1}{2}(\theta - \theta_0)^{\top} H(\theta_0)(\theta - \theta_0)$ describes the curvature effects. The eigenvalues of $H(\theta_0)$ dictate the
nature of the critical point $\theta_0$. Specifically, if all $\lambda_i > 0$, $\theta_0$ is a local minimum; if all $\lambda_i < 0$, it is a local maximum; and if the eigenvalues have mixed signs, $\theta_0$ is a saddle point. The leading-order approximation to the change in the loss function, $\Delta\mathcal{L} \approx \frac{1}{2}\delta\theta^\top H(\theta_0)\delta\theta$, highlights the dependence on the eigenstructure of $H(\theta_0)$. In the context of gradient descent, parameter updates follow the iterative scheme:

$$ \theta^{(t+1)} = \theta^{(t)} - \eta \nabla_{\theta} \mathcal{L}(\theta^{(t)}), \qquad (3.21) $$

where $\eta$ is the learning rate. Substituting the Taylor expansion of $\nabla_{\theta} \mathcal{L}(\theta^{(t)})$ around $\theta_0$ gives:

$$ \theta^{(t+1)} = \theta^{(t)} - \eta \left[ \nabla_{\theta} \mathcal{L}(\theta_0) + H(\theta_0)(\theta^{(t)} - \theta_0) \right]. \qquad (3.22) $$

To analyze this update rigorously, we project $\theta^{(t)} - \theta_0$ onto the eigenbasis of $H(\theta_0)$, expressing it as:

$$ \theta^{(t)} - \theta_0 = \sum_{i=1}^{d} c_i^{(t)} \mathbf{v}_i, \qquad (3.23) $$

where $c_i^{(t)} = \mathbf{v}_i^\top (\theta^{(t)} - \theta_0)$. Substituting this expansion into the gradient descent update rule yields:

$$ c_i^{(t+1)} = c_i^{(t)} - \eta \left[ \mathbf{v}_i^\top \nabla_{\theta} \mathcal{L}(\theta_0) + \lambda_i c_i^{(t)} \right]. \qquad (3.24) $$

The convergence of this iterative scheme is governed by the condition $|1 - \eta\lambda_i| < 1$, which constrains the learning rate $\eta$ relative to the spectrum of $H(\theta_0)$. For eigenvalues $\lambda_i$ with large magnitudes, excessively large learning rates $\eta$ can cause oscillatory or divergent updates. In the Neural Tangent

<figure><img src="image_33.png" /><figcaption>Figure 3.3: Hessian Structure (NTK Regime) & NTK Matrix</figcaption></figure>

Kernel (NTK) regime, the evolution of a neural network during training can be approximated by a linearization of the network output around the initialization. Let $f_\theta(x)$ denote the output of the network for input $x$. Linearizing $f_\theta(x)$ around $\theta_0$ gives:

$$ f_{\theta}(x) \approx f_{\theta_0}(x) + \nabla_{\theta} f_{\theta_0}(x)^{\top} (\theta - \theta_0). \qquad (3.25) $$

The NTK, defined as:

$$ K(x, x') = \nabla_{\theta} f_{\theta_0}(x)^{\top} \nabla_{\theta} f_{\theta_0}(x'), \qquad (3.26) $$

remains approximately constant during training for sufficiently wide networks. The training dynamics of the parameters are described by:

$$ \frac{d\theta}{dt} = -\nabla_{\theta} \mathcal{L}(\theta), \qquad (3.27) $$
<figure><img src="image_34.png" /><figcaption>Figure 3.4: Hessian Spectrum (NTK Regime)</figcaption></figure>

which, under the NTK approximation, becomes:

$$ \frac{d\theta}{dt} = -K\nabla_{\theta}\mathcal{L}(\theta), \qquad (3.28) $$

where $K$ is the NTK matrix evaluated at initialization. The evolution of the loss function is governed by the eigenvalues of $K$, which control the rate of convergence in different directions.

The spectral properties of the Hessian play a pivotal role in the generalization properties of neural networks. Empirical studies reveal that the eigenvalue spectrum of $H(\theta)$ often exhibits a "bulk-and-spike" structure, with a dense bulk of eigenvalues near zero and a few large outliers. The bulk corresponds to flat directions in the loss landscape, which contribute to the robustness and generalization of the model, while the spikes represent sharp directions associated with overfitting. This spectral structure can be analyzed using random matrix theory, where the density of eigenvalues $\rho(\lambda)$ is modeled by distributions such as the Marchenko-Pastur law:

$$ \rho(\lambda) = \frac{1}{2\pi\lambda q} \sqrt{(\lambda_+ - \lambda)(\lambda - \lambda_-)}, \qquad (3.29) $$

where $\lambda_{\pm} = (1 \pm \sqrt{q})^2$ are the spectral bounds and $q = \frac{d}{n}$ is the ratio of the number of parameters to the number of data points. This rigorous analysis links the Hessian structure to both the optimization dynamics and the generalization performance of neural networks, providing a comprehensive mathematical understanding of the training process. The Hessian $H(\theta)$ satisfies:

$$ H(\theta) = \nabla_{\theta}^2 \mathcal{L}(\theta) = \mathbb{E}_{(x,y)} \left[ \nabla_{\theta} f_{\theta}(x) \nabla_{\theta} f_{\theta}(x)^{\top} \right]. \qquad (3.30) $$

For overparameterized networks, $H(\theta)$ is nearly degenerate, implying the existence of flat minima.
### 3.3.4 NTK Linearization

The dynamics of neural networks under gradient flow can be comprehensively described by beginning with the parameterized representation of the network $f_\theta(x)$, where $\theta \in \mathbb{R}^p$ denotes the set of trainable parameters, $x \in \mathbb{R}^d$ is the input, and $f_\theta(x) \in \mathbb{R}^m$ represents the output. The objective of training is to minimize a loss function $\mathcal{L}(\theta)$, defined over a dataset $\{(x_i, y_i)\}_{i=1}^n$, where $x_i \in \mathbb{R}^d$ and $y_i \in \mathbb{R}^m$ represent the input-target pairs. The evolution of the parameters during training is governed by the gradient flow equation $\frac{d\theta}{dt} = -\nabla_\theta \mathcal{L}(\theta)$, where $\nabla_\theta \mathcal{L}(\theta)$ is the gradient of the loss function with respect to the parameters. To analyze the dynamics of the network outputs, we first consider the time derivative of $f_\theta(x)$. Using the chain rule, this is expressed as:

$$ \frac{\partial f_{\theta}(x)}{\partial t} = \nabla_{\theta} f_{\theta}(x)^{\top} \frac{d\theta}{dt}. \qquad (3.31) $$

Substituting $\frac{d\theta}{dt} = -\nabla_{\theta} \mathcal{L}(\theta)$, we have:

$$ \frac{\partial f_{\theta}(x)}{\partial t} = -\nabla_{\theta} f_{\theta}(x)^{\top} \nabla_{\theta} \mathcal{L}(\theta). \qquad (3.32) $$

The gradient of the loss function, $\mathcal{L}(\theta)$, can be expressed explicitly in terms of the training data. For a generic loss function over the dataset, this takes the form:

$$ \mathcal{L}(\theta) = \frac{1}{n} \sum_{i=1}^{n} \ell(f_{\theta}(x_i), y_i), \qquad (3.33) $$

where $\ell(f_\theta(x_i), y_i)$ represents the loss for the $i$-th data point. The gradient of the loss with respect to the parameters is therefore given by:

$$ \nabla_{\theta} \mathcal{L}(\theta) = \frac{1}{n} \sum_{i=1}^{n} \nabla_{\theta} f_{\theta}(x_i) \nabla_{f_{\theta}(x_i)} \ell(f_{\theta}(x_i), y_i). \qquad (3.34) $$

Substituting this back into the time derivative of $f_\theta(x)$, we obtain:

$$ \frac{\partial f_{\theta}(x)}{\partial t} = -\frac{1}{n} \sum_{i=1}^{n} \nabla_{\theta} f_{\theta}(x)^{\top} \nabla_{\theta} f_{\theta}(x_i) \nabla_{f_{\theta}(x_i)} \ell(f_{\theta}(x_i), y_i). \qquad (3.35) $$

To introduce the Neural Tangent Kernel (NTK), we define it as the Gram matrix of the Jacobians of the network output with respect to the parameters:

$$ \Theta(x, x'; \theta) = \nabla_{\theta} f_{\theta}(x)^{\top} \nabla_{\theta} f_{\theta}(x'). \qquad (3.36) $$

Using this definition, the time evolution of the output becomes:

$$ \frac{\partial f_{\theta}(x)}{\partial t} = -\frac{1}{n} \sum_{i=1}^{n} \Theta(x, x_i; \theta) \nabla_{f_{\theta}(x_i)} \ell(f_{\theta}(x_i), y_i). \qquad (3.37) $$

In the overparameterized regime, where the number of parameters $p$ is significantly larger than the number of training data points $n$, it has been empirically and theoretically observed that the NTK $\Theta(x, x'; \theta)$ remains nearly constant during training. Specifically, $\Theta(x, x'; \theta) \approx \Theta(x, x'; \theta_0)$, where $\theta_0$ represents the parameters at initialization. This constancy significantly simplifies the analysis of the network's training dynamics. To see this, consider the solution to the differential equation governing the output dynamics. Let $F(t) \in \mathbb{R}^{n \times m}$ represent the matrix of network outputs for all training inputs, where the $i$-th row corresponds to $f_\theta(x_i)$. The dynamics can be expressed in matrix form as:

$$ \frac{\partial F(t)}{\partial t} = -\frac{1}{n} \Theta(\theta_0) \nabla_{F} \mathcal{L}(F), \qquad (3.38) $$
<figure><img src="image_35.png" /><figcaption>Figure 3.5: Final Fit & Training Dynamics (NN vs NTK)</figcaption></figure>

where $\Theta(\theta_0) \in \mathbb{R}^{n \times n}$ is the NTK matrix evaluated at initialization, and $\nabla_F \mathcal{L}(F)$ is the gradient of the loss with respect to the output matrix $F$. For the special case of a mean squared error loss, $\mathcal{L}(F) = \frac{1}{2n}\|F - Y\|_F^2$, where $Y \in \mathbb{R}^{n \times m}$ is the matrix of target outputs, the gradient simplifies to:

$$ \nabla_F \mathcal{L}(F) = \frac{1}{n}(F - Y). \qquad (3.39) $$

Substituting this into the dynamics, we obtain:

$$ \frac{\partial F(t)}{\partial t} = -\frac{1}{n^2}\Theta(\theta_0)(F(t) - Y). \qquad (3.40) $$

The solution to this differential equation is:

$$ F(t) = Y + e^{-\frac{\Theta(\theta_0)}{n^2}t}(F(0) - Y), \qquad (3.41) $$

where $F(0)$ represents the initial outputs of the network. As $t \to \infty$, the exponential term vanishes, and the network outputs converge to the targets $Y$, provided that $\Theta(\theta_0)$ is positive definite. The rate of convergence is determined by the eigenvalues of $\Theta(\theta_0)$, with smaller eigenvalues corresponding to slower convergence along the associated eigenvectors. To understand the stationary points of this system, we note that these occur when $\frac{\partial F(t)}{\partial t} = 0$. From the dynamics, this implies:

$$ \Theta(\theta_0)(F - Y) = 0. \qquad (3.42) $$

If $\Theta(\theta_0)$ is invertible, this yields $F = Y$, indicating that the network exactly interpolates the training data at the stationary point. However, if $\Theta(\theta_0)$ is not full-rank, the stationary points form a subspace of solutions satisfying $(I - \Pi)(F - Y) = 0$, where $\Pi$ is the projection operator onto the column space of $\Theta(\theta_0)$.

The NTK framework provides a mathematically rigorous lens to analyze training dynamics, elucidating the interplay between parameter evolution, kernel properties, and loss convergence in neural networks. By linearizing the training dynamics through the NTK, we achieve a deep understanding of how overparameterized networks evolve under gradient flow and how they reach stationary points, revealing their capacity to interpolate data with remarkable precision.
3.3.5 Python Code to generate Figure 3.2 illustrating the Gradient Flow and Stationary Points

The Python code below produces the Figure 3.2 illustrating the Gradient Flow and Stationary Points.

import numpy as np
import matplotlib.pyplot as plt
import sympy as sp

# ==============================
# 1. Define the function f(x,y)
# ==============================
x, y = sp symbols('x y')
f = x**3 - 3*x*y**2 # Example: a saddle-like function (Monkey Saddle)

# Gradient
grad_f = [sp.diff(f, var) for var in (x, y)]

# ==============================
# 2. Solve for stationary points
# ==============================
stationary_points = sp solve Object
stationary_points = [(float(p[x]), float(p[y])) for p in
                       stationary_points]

print("Stationary Points:", stationary_points)

# ==============================
# 3. Convert symbolic to numerical fn
# ==============================
f_num = sp lambdify((x, y), f, 'numpy')
grad_f_num = [sp lambdify((x, y), grad, 'numpy') for grad in grad_f]

# ==============================
# 4. Create meshgrid for vector field
# ==============================
X, Y = np meshgrid(np linspace(-2, 2, 30), np linspace(-2, 2, 30))
U = grad_f_num[0](X, Y) # df/dx
V = grad_f_num[1](X, Y) # df/dy

# Normalize arrows for better visualization
N = np sqrt(U**2 + V**2)
U, V = U/N, V/N

# ==============================
# 5. Plot gradient flow
# ==============================
plt figure (figsize=(7,6))

# Function contour
plt contour(X, Y, f_num(X, Y), levels=20, cmap="viridis")

# Gradient flow (vector field)
plt quiver(X, Y, -U, -V, color="red", alpha=0.6) # negative for
```
gradient descent flow

# Stationary points
for (sx, sy) in stationary_points:
    plt.plot(sx, sy, 'bo', markersize=10, label="Stationary Point")

plt.title("Gradient Flow and Stationary Points")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()

plt.grid(True)
plt.show()
```

### 3.3.6 Python Code to generate Figures 3.3 and 3.4 illustrating the Hessian Structure

The Python code below produces the Figures 3.3 and 3.4 illustrating the Hessian Structure.

```python
import torch
import torch.nn as nn
import matplotlib.pyplot as plt

# =============================
# 1. Define a simple MLP
# =============================
class Net nn.Module):
    def __init__(self, d_in, d_hidden, d_out):
        super().__init__()
        self.fc1 = nn.Linear(d_in, d_hidden)
        self.fc2 = nn.Linear(d_hidden, d_out)

    def forward(self, x):
        return torch.tanh(self.fc1(x)) @ self.fc2.weight.T + self.fc2.bias

# =============================
# 2. Compute NTK
# =============================
def compute_ntk(model, x):
    f = model(x)
    grads = []

    for i in range(f.shape[0]):
        grad = torch.autograd.grad(f[i, 0], model.parameters(),
                                 retain_graph=True, create_graph=True)
        grads.append torch.cat([g Reshape(-1) for g in grad]))

    grads = torch.stack grads)
    NTK = grads @ grads.T
    return NTK Detach()

# =============================
# 3. Compute Hessian
# =============================
```
```python
def compute_hessian(loss, params):
    grads = torch.autograd.grad(loss, params, create_graph=True)
    grads_flat = torch.cat([g.reshape(-1) for g in grads])
    n = grads_flat.shape[0]
    H = torch(zeros(n, n))
    for i in range(n):
        grad2 = torch.autograd.grad grads_flat[i], params, retain_graph=True)
        grad2_flat = torch.cat([g.reshape(-1) for g in grad2])
        H[i] = grad2_flat
    return H

# =============================
# 4. Dataset
# =============================
    torch.manual_seed(0)
    N = 30
    x = torch.linspace(-3, 3, N).reshape(-1,1)
    y_true = torch.sin(x)

# =============================
# 5. Initialize wide network (NTK regime)
# =============================
    net = Net(1, 200, 1)
    params = list(net.parameters())

# Forward pass and loss
    y_pred = net(x)
    loss = nn.MSELoss() (y_pred, y_true)

# =============================
# 6. Compute NTK and Hessian
# =============================
    NTK = compute_ntk(net, x).numpy()
    H = compute_hessian(loss, params).detach(). numpy()

# =============================
# 7. Plot Hessian and NTK side by side
# =============================
    plt.figure(figsize=(12,5))
    plt.subplots(1,2,1)
    plt.imshow(H, cmap='bwr')
    plt.colorbar()
    plt.title('Hessian Structure (NTK Regime)')
    plt>xlabel('Parameter Index')
    plt>ylabel('Parameter Index')

# =============================
# 8. Plot NTK Matrix
# =============================
    plt.subplots(1,2,1)
    plt.imshow(NTK, cmap='bwr')
    plt.colorbar()
    plt.title('NTK Matrix')
    plt>xlabel('Data Point Index')
    plt>ylabel('Data Point Index')
```
```python
89 plt.tight_layout()
90 plt.show()
91
92 # ===================
93 # 8. Optional: Hessian spectrum
94 # ===================
95 eigvals_H = torch.laloog.eigvals(torch.tensor(H)).real
96 pltsemiology(sorted(eigvals_H, reverse=True))
97 plt.title("Hessian Spectrum (NTK Regime)")
98 plt.xlabel("Eigenvalue Index")
99 plt.ylabel("Eigenvalue")
100 plt.show()
```

### 3.3.7 Python Code to generate Figure 3.5 illustrating the Final Fit & Training Dynamics (NN vs NTK)

The Python code below produces the Figure 3.5 illustrating the Final Fit & Training Dynamics (NN vs NTK).

```python
1 import torch
2 import torch.nn as nn
3 import torch.nn(functional as F
4 import matplotlib.pyplot as plt
5
6 # ===================
7 # 1. Define a simple MLP
8 # ===================
9 class Net(nn.Module):
10     def __init__(self, d_in, d_hidden, d_out):
11         super().__init__()
12         self.fc1 = nn.Linear(d_in, d_hidden)
13         self.fc2 = nn.Linear(d_hidden, d_out)
14
15     def forward(self, x):
16         x = torch.tanh(self.fc1(x))
17         x = self.fc2(x)
18         return x
19
20 # ===================
21 # 2. Compute NTK
22 # ===================
23 def compute_ntk(model, x):
24     """
25     Compute NTK matrix
26     """
27     f = model(x)
28     grads = []
29
30     for i in range(f.shape[0]):
31         grad = torch.autograd.grad(f[i, 0], model.parameters(),
32                                 retain_graph=True, create_graph=True)
33         grads.append MissCat([greshape(-1) for g in grad]))
34     grads = torch.stack grads) # [N, num_params]
35```
```python
35     NTK = grads @ grads.T          # Gram matrix
36     return NTK detach()
37
38     # ======================
39     # 3. Train and compare
40     # ======================
41   torch.manual_seed(0)
42   N = 20
43   x = torch.linspace(-3, 3, N).reshape(-1, 1)
44   y_true = torch.sin(x)          # target function
45
46   net = Net(1, 50, 1)
47   loss_fn = nn.MSELoss()
48   opt = torch.optim.SGD(net.parameters(), lr=0.05)
49
50     # NTK at initialization
51   NTK = compute_ntk(net, x)
52   f0 = net(x).detach()
53   theta_inv = torch.lerpg.nv(NTK)   # pseudo-inverse
54   alpha = theta_inv @ (y_true - f0)  # coefficients for NTK linearization
55   f_ntk = f0.clone()
56
57     # Store training trajectories
58   preds_nn, preds_ntk = [], []
59
60   for step in range(200):
61     # --- Full NN training ---
62     opt.zeroGrad()
63     y_pred = net(x)
64     loss = loss_fn(y_pred, y_true)
65     loss backward()
66     opt.step()
67
68     preds_nn.append(y_pred Detach().clone())
69
70     # --- NTK linearized dynamics (gradient flow approx) ---
71     decay = torch.exp IDetach().float32())
72     f_ntk = f0 + NTK @ alpha * (1 - decay)
73     preds_ntk.append(f_ntk.clone())
74
75     # ======================
76     # 4. Plot results
77     # ======================
78   plt.figure(figsize=(12, 5))
79
80     # Final fits
81   plt.subplots(1, 2, 1)
82   plt.scatter x, y_true, label="True Function", c="black")
83   plt.plot x, preds_nn[-1], label="NN trained", c="blue")
84   plt.plot x, preds_ntk[-1], label="NTK linearization", c="red", linewidth=2)
85
86   plt.legend()
87   plt.title("Final Fit")
88
```
```python
# Training trajectories (first 50 steps)
plt.subplots(1,2,2)
for t in [0, 5, 10, 20, 50, 100, 199]:
    plt.plot(x, preds_nn[t], 'o', markersize=4, alpha=0.5, label=f"NN step {t}")
    if t in [0,199] else None, c="blue")
    plt.plot(x, preds_ntk[t], 'o', markersize=4, alpha=0.5, label=f"NTK step {t}")
    if t in [0,199] else None, c="red")

plt.scatter(x, y_true, c="black", s=20)
plt.legend()
plt.title("Training Dynamics (NN vs NTK)")
plt.show()
```

## 3.4 NTK Regime

### 3.4.1 Literature Review of NTK Regime

Jacot et. al. (2018) [94] in a seminal paper introduced the Neural Tangent Kernel (NTK) and establishes its theoretical foundation. The authors show that in the infinite-width limit, the dynamics of gradient descent in neural networks can be described by a kernel method, where the NTK remains constant during training. This work bridges the gap between deep learning and kernel methods, providing a framework to analyze the training and generalization of wide neural networks. Lee et. al. (2017) [95] did a work that predates the NTK but lays the groundwork by showing that infinitely wide neural networks behave as Gaussian processes. The authors derive the kernel corresponding to such networks, which is a precursor to the NTK. This paper is crucial for understanding the connection between neural networks and kernel methods. Chizat and Bach (2018) [484] provided a rigorous analysis of gradient descent in over-parameterized models, including neural networks. It complements the NTK framework by showing that gradient descent converges to global minima in such settings. The work highlights the role of over-parameterization in simplifying the optimization landscape. Du et. al. (2019) [485] proved that gradient descent can find global minima in deep neural networks under the NTK regime. The authors provide explicit convergence rates and show that the NTK framework guarantees efficient optimization for wide networks. This work strengthens the theoretical understanding of why deep learning works. Arora et. al. (2019) [486] provided a fine-grained analysis of optimization and generalization in two-layer neural networks under the NTK regime. It establishes precise bounds on the generalization error and shows how over-parameterization leads to benign optimization landscapes. Allen-Zhu et. al. (2019) [487] extended the NTK framework to deep networks and provided a comprehensive convergence theory for over-parameterized neural networks. The authors show that gradient descent converges to global minima and that the NTK remains approximately constant during training. Cao and Gu (2019) [488] derived generalization bounds for wide and deep neural networks trained with stochastic gradient descent (SGD) under the NTK regime. It highlights the role of the NTK in controlling the generalization error and provides insights into the implicit regularization of SGD. Yang (2019) [489] generalized the NTK framework to architectures with weight sharing, such as convolutional neural networks (CNNs). The author derives the NTK for such architectures and shows that they also exhibit Gaussian process behavior in the infinite-width limit. Huang and Yau (2020) [490] extended the NTK framework by introducing the Neural Tangent Hierarchy (NTH), which captures higher-order interactions in the training dynamics of deep networks. The authors provide a more refined analysis of the training process beyond the first-order approximation of the NTK. Belkin et. al. (2019) [491] explored the connection between deep learning and kernel learning, emphasizing the role of the NTK in understanding generalization and optimization. It
provides a high-level perspective on why the NTK framework is essential for analyzing modern machine learning models.

### 3.4.2 Analysis of NTK Regime

The Neural Tangent Kernel (NTK) regime is a fundamental framework for understanding the dynamics of gradient descent in highly overparameterized neural networks. Consider a neural network $f(\mathbf{x}; \boldsymbol{\theta})$ parameterized by $\boldsymbol{\theta} \in \mathbb{R}^P$, where $P$ represents the total number of parameters, and $\mathbf{x} \in \mathbb{R}^d$ is the input vector. For a training dataset $\{(\mathbf{x}_i, y_i)\}_{i=1}^N$, the loss function $L(t)$ at time $t$ is given by

$$L(t) = \frac{1}{2N} \sum_{i=1}^{N} (f(\mathbf{x}_i; \boldsymbol{\theta}(t)) - y_i)^2. \qquad (3.43)$$

The parameters evolve according to gradient descent as $\boldsymbol{\theta}(t+1) = \boldsymbol{\theta}(t) - \eta \nabla_{\boldsymbol{\theta}} L(t)$, where $\eta > 0$ is the learning rate. In the NTK regime, we consider the first-order Taylor expansion of the network

<figure><img src="image_36.png" /><figcaption>Figure 3.6: Final Predictions in NTK Regime & Training Dynamics (NN vs NTK)</figcaption></figure>

output around the initialization $\boldsymbol{\theta}_0$:

$$f(\mathbf{x}; \boldsymbol{\theta}) \approx f(\mathbf{x}; \boldsymbol{\theta}_0) + \nabla_{\boldsymbol{\theta}} f(\mathbf{x}; \boldsymbol{\theta}_0)^{\top} (\boldsymbol{\theta} - \boldsymbol{\theta}_0). \qquad (3.44)$$

This linear approximation transforms the nonlinear dynamics of $f$ into a simpler, linearized form. To analyze training, we introduce the Jacobian matrix $J \in \mathbb{R}^{N \times P}$, where $J_{ij} = \frac{\partial f(\mathbf{x}_i; \boldsymbol{\theta}_0)}{\partial \boldsymbol{\theta}_j}$. The vector of outputs $\mathbf{f}(t) \in \mathbb{R}^N$, aggregating predictions over the dataset, evolves as

$$\mathbf{f}(t) = \mathbf{f}(0) + J(\boldsymbol{\theta}(t) - \boldsymbol{\theta}_0). \qquad (3.45)$$

The NTK $\Theta \in \mathbb{R}^{N \times N}$ is defined as

$$\Theta_{ij} = \nabla_{\boldsymbol{\theta}} f(\mathbf{x}_i; \boldsymbol{\theta}_0)^{\top} \nabla_{\boldsymbol{\theta}} f(\mathbf{x}_j; \boldsymbol{\theta}_0). \qquad (3.46)$$

As $P \to \infty$, the NTK converges to a deterministic matrix that remains nearly constant during training. Substituting the linearized form of $\mathbf{f}(t)$ into the gradient descent update equation gives

$$\mathbf{f}(t+1) = \mathbf{f}(t) - \frac{\eta}{N} \Theta(\mathbf{f}(t) - \mathbf{y}), \qquad (3.47)$$
where $\mathbf{y} \in \mathbb{R}^N$ is the vector of true labels. Defining the residual $\mathbf{r}(t) = \mathbf{f}(t) - \mathbf{y}$, the dynamics of training reduce to

$$ \mathbf{r}(t+1) = \left(I - \frac{\eta}{N}\Theta\right) \mathbf{r}(t). \qquad (3.48) $$

The eigendecomposition $\Theta = Q\Lambda Q^\top$, with orthogonal $Q$ and diagonal $\Lambda = \text{diag}(\lambda_1, \dots, \lambda_N)$, allows us to analyze the decay of residuals in the eigenbasis of $\Theta$:

$$ \tilde{\mathbf{r}}(t+1) = \left(I - \frac{\eta}{N}\Lambda\right) \tilde{\mathbf{r}}(t), \qquad (3.49) $$

where $\tilde{\mathbf{r}}(t) = Q^\top \mathbf{r}(t)$. Each component decays as

$$ \tilde{r}_i(t) = \left(1 - \frac{\eta\lambda_i}{N}\right)^t \tilde{r}_i(0). \qquad (3.50) $$

For small $\eta$, the training dynamics are approximately continuous, governed by

$$ \frac{d\mathbf{r}(t)}{dt} = -\frac{1}{N}\Theta\mathbf{r}(t), \qquad (3.51) $$

leading to the solution

$$ \mathbf{r}(t) = \exp\left(-\frac{\Theta t}{N}\right) \mathbf{r}(0). \qquad (3.52) $$

The NTK for specific architectures, such as fully connected ReLU networks, can be derived using layerwise covariance matrices. Let $\Sigma^{(l)}(\mathbf{x}, \mathbf{x}')$ denote the covariance between pre-activations at layer $l$. The recurrence relation for $\Sigma^{(l)}$ is

$$ \Sigma^{(l)}(\mathbf{x}, \mathbf{x}') = \frac{1}{2\pi} \|\mathbf{z}^{(l-1)}(\mathbf{x})\| \|\mathbf{z}^{(l-1)}(\mathbf{x}')\| (\sin \theta + (\pi - \theta) \cos \theta), \qquad (3.53) $$

where $\theta = \cos^{-1}\left(\frac{\Sigma^{(l-1)}(\mathbf{x}, \mathbf{x}')}{\sqrt{\Sigma^{(l-1)}(\mathbf{x}, \mathbf{x})\Sigma^{(l-1)}(\mathbf{x}', \mathbf{x}')}}\right)$. The NTK, a sum over contributions from all layers, quantifies how parameter updates propagate through the network.

In the infinite-width limit, the NTK framework predicts generalization properties, as the kernel matrix $\Theta$ governs both training and test-time behavior. The NTK connects neural networks to classical kernel methods, offering a bridge between deep learning and well-established theoretical tools in approximation theory. This regime's deterministic and analytical tractability enables precise characterizations of network performance, convergence rates, and robustness to initialization and learning rate variations.

### 3.4.3 Python Code to generate Figure 3.6 illustrating the Final Predictions in NTK Regime & Training Dynamics (NN vs NTK)

The Python code below produces the Figure 3.6 illustrating the Final Predictions in NTK Regime & Training Dynamics (NN vs NTK).

```python
import torch
import torch.nn as nn
import matplotlib.pyplot as plt

# ===================
# 1. Simple MLP
# ===================
class Net nn.Module):
    def __init__(self, d_in, d_hidden, d_out):
```
```python
super().__init__()
self.fc1 = nn.Linear(d_in, d_hidden)
self.fc2 = nn.Linear(d_hidden, d_out)

def forward(self, x):
    return torch.tanh(self.fc1(x)) @ self.fc2.weight.T + self.fc2.bias  # simple tanh MLP

# ==============================
# 2. Compute NTK at initialization
# ==============================
def compute_ntk(model, x):
    f = model(x)
    grads = []

    for i in range(f.shape[0]):
        grad = torch.autograd.grad(f[i, 0], model.parameters(),
                                 retain_graph=True, create_graph=True)
        grads.append MissedCat([greshape(-1) for g in grad]))

    grads = torch.stack grads)
    NTK = grads @ grads.T
    return NTK Detach()

# ==============================
# 3. Prepare dataset
# ==============================
torch.manual_seed(0)
N = 50
x = torch.linspace(-3, 3, N).reshape(-1, 1)
y_true = torch.sin(x)

# ==============================
# 4. Initialize model
# ==============================
net = Net(1, 200, 1) # large hidden layer to NTK regime
loss_fn = nn.MSELoss()
opt = torch.optim.SGD(net.parameters(), lr=0.05)

# ==============================
# 5. NTK linearization
# ==============================
f0 = net(x).detach()
NTK = compute_ntk(net, x)
theta_inv = torch.linalg.pnv(NTK)
alpha = theta_inv @ (y_true - f0)

# ==============================
# 6. Training and NTK evolution
# ==============================
steps = 200
predsnn, predsntk = [], []

for step in range(steps):
    # Full NN training
```
```python
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
```
# 4 Generalization Bounds: PAC-Bayes and Spectral Analysis

The PAC-Bayes framework provides a rigorous mathematical approach to deriving generalization bounds for machine learning models, particularly in stochastic settings. Given a hypothesis space $\mathcal{H}$, a prior distribution $P$ over $\mathcal{H}$ (encoding initial beliefs before observing data), and a posterior distribution $Q$ (learned from data), the PAC-Bayes theorem bounds the expected generalization error. Let $\mathcal{L}(h)$ denote the true risk of hypothesis $h \in \mathcal{H}$ and $\hat{\mathcal{L}}(h)$ its empirical risk on a sample $S$ of size $m$. The key PAC-Bayes bound states that, with probability at least $1 - \delta$ over the draw of $S$, the following holds for all posteriors $Q$:

$$ \mathbb{E}_{h \sim Q}[\mathcal{L}(h)] \le \mathbb{E}_{h \sim Q}[\hat{\mathcal{L}}(h)] + \sqrt{\frac{KL(Q \| P) + \log \frac{m}{\delta}}{2(m-1)}}, \quad (4.1) $$

where $KL(Q \| P)$ is the Kullback-Leibler divergence between $Q$ and $P$. This bound balances empirical performance and the complexity of the posterior relative to the prior, ensuring generalization when the learned distribution $Q$ does not deviate excessively from $P$. The theorem is particularly useful in Bayesian deep learning and stochastic optimization, where $Q$ is often taken as a Gaussian perturbation of a deterministic model.

Spectral analysis examines generalization through the lens of the loss landscape and the eigen-spectrum of the Hessian or the Neural Tangent Kernel (NTK). Let $\mathbf{H}$ be the Hessian of the loss at a local minimum, and let $\{\lambda_i\}$ denote its eigenvalues. The effective dimensionality of the solution, which influences generalization, can be characterized by the ratio of the trace to the largest eigenvalue:

$$ R = \frac{\mathrm{Tr}(\mathbf{H})}{\lambda_{\max}}. \quad (4.2) $$

Sharp minima, associated with poor generalization, correspond to large $\lambda_{\max}$ relative to the trace, while flat minima (with smaller $\lambda_{\max}$) generalize better. Additionally, the NTK's eigenspectrum determines the speed of gradient descent convergence and the stability of learning dynamics. Let $\mathbf{H}_0$ be the NTK at initialization; its eigenvalue decay rate governs how quickly different modes of the target function are learned. If $\mathbf{H}_0$ has a rapidly decaying spectrum, high-frequency components (corresponding to small eigenvalues) are learned slowly, leading to implicit regularization. The generalization error can be bounded in terms of the NTK's eigendecomposition, linking spectral properties to model performance. This approach provides a mechanistic understanding of why overparameterized models generalize despite their capacity to overfit.

## 4.1 PAC-Bayes Formalism

### 4.1.1 Literature Review of PAC-Bayes Formalism

McAllester (1999) [99] introduced the PAC-Bayes bound, a fundamental theorem that provides generalization guarantees for Bayesian learning models. He established a trade-off between com-
plexity and empirical risk, serving as the theoretical foundation for modern PAC-Bayesian analysis. Catoni (2007) [100] in his book rigorously extended the PAC-Bayes framework by linking it with information-theoretic and statistical mechanics concepts and introduced exponential and Gibbs priors for learning, improving PAC-Bayesian bounds for supervised classification. Germain et. al. (2009) [101] applied PAC-Bayes theory to linear classifiers, including SVMs and logistic regression. They demonstrated that PAC-Bayesian generalization bounds are tighter than classical Vapnik-Chervonenkis (VC) dimension bounds. Seeger (2002) [102] extended PAC-Bayes bounds to Gaussian Process models, proving tight generalization guarantees for Bayesian classifiers. He laid the groundwork for probabilistic kernel methods. Alquier et. al. (2006) [103] connected variational inference and PAC-Bayes bounds, proving that variational approximations can preserve generalization guarantees of PAC-Bayesian bounds. Dziugaite and Roy (2017) [104] gave one of the first applications of PAC-Bayes to deep learning. They derived nonvacuous generalization bounds for stochastic neural networks, bridging theory and practice. Cox and Ghosh (2022) [54] introduce an analogue of the Mertens function that exhibits similar oscillatory behavior, analyze its properties using Shapiro's theorem, generalized Dirichlet products, and generalized Möbius inversion, thereby offering new insights into the analytic structure underlying Mertens-type functions. Rivasplata et. al. (2020) [106] provided novel PAC-Bayes bounds that improve over existing guarantees, making PAC-Bayesian bounds more practical for modern ML applications. Lever et. al. (2013) [107] explored data-dependent priors in PAC-Bayes theory, showing that adaptive priors lead to tighter generalization bounds. Rivasplata et. al. (2018) [108] introduced instance-dependent priors, improving personalized learning and making PAC-Bayesian methods more useful for real-world machine learning problems. Lindemann et. al. (2024) [109] integrated PAC-Bayes theory with conformal prediction to improve formal verification in control systems, demonstrating PAC-Bayes' relevance to safety-critical applications.

## 4.1.2 Analysis of PAC-Bayes Formalism

The PAC-Bayes formalism is a foundational framework in statistical learning theory, designed to provide probabilistic guarantees on the generalization performance of learning algorithms. By combining principles from the PAC (Probably Approximately Correct) framework and Bayesian reasoning, PAC-Bayes delivers bounds that characterize the expected performance of hypotheses drawn from posterior distributions, given a finite sample of data. This document presents an extremely rigorous and mathematically precise description of the PAC-Bayes formalism, emphasizing its theoretical constructs and implications. At the core of the PAC-Bayes formalism lies the ambition to rigorously quantify the generalization ability of hypotheses $h \in \mathcal{H}$ based on their performance on a finite dataset $S \sim \mathcal{D}^m$, where $\mathcal{D}$ represents the underlying, and typically unknown, data distribution. The PAC framework, which was originally designed to provide high-confidence guarantees on the true risk

$$R(h) = \mathbb{E}_{z \sim \mathcal{D}}[\ell(h, z)], \quad (4.3)$$

is enriched in PAC-Bayes by incorporating principles from Bayesian reasoning. This integration allows for bounds not just on individual hypotheses but on distributions $Q$ over $\mathcal{H}$, yielding a sophisticated characterization of generalization that inherently accounts for the variability and uncertainty in the hypothesis space. There are some Mathematical Constructs: True and Empirical Risks. The true risk $R(h)$, as defined by the expected loss, is typically inaccessible due to the unknown nature of $\mathcal{D}$. Instead, the empirical risk

$$\hat{R}(h, S) = \frac{1}{m} \sum_{i=1}^{m} \ell(h, z_i) \quad (4.4)$$

serves as a computable proxy. The key question addressed by PAC-Bayes is: *How does $\hat{R}(h, S)$ relate to $R(h)$, and how can we bound the deviation probabilistically?* For a distribution $Q$ over $\mathcal{H}$,
<figure><img src="image_37.png" /><figcaption>Figure 4.1: PAC-Bayes Bound vs. Sample Size (KL(Q||P) = 2.0)</figcaption></figure>

these risks are generalized as:

$$R(Q) = \mathbb{E}_{h \sim Q}[R(h)], \quad \hat{R}(Q, S) = \mathbb{E}_{h \sim Q}[\hat{R}(h, S)]. \qquad (4.5)$$

This generalization is pivotal because it allows the analysis to transcend individual hypotheses and consider probabilistic ensembles, where $Q(h)$ represents a posterior belief over the hypothesis space conditioned on the observed data. We now need to discuss how Prior and Posterior Distributions encode knowledge and complexity. The prior $P$ is a fixed distribution over $\mathcal{H}$ that reflects pre-data assumptions about the plausibility of hypotheses. Crucially, $P$ must be independent of $S$ to avoid biasing the bounds. The posterior $Q$, however, is data-dependent and typically chosen to minimize a combination of empirical risk and complexity. This choice is guided by the PAC-Bayes inequality, which regularizes $Q$ via its Kullback-Leibler (KL) divergence from $P$:

$$\mathrm{KL}(Q \| P) = \int_{\mathcal{H}} Q(h) \log \frac{Q(h)}{P(h)} dh. \qquad (4.6)$$

The KL divergence quantifies the informational cost of updating $P$ to $Q$, serving as a penalty term that discourages overly complex posteriors. This regularization is critical in preventing overfitting, ensuring that $Q$ achieves a balance between data fidelity and model simplicity.

Let's derive the PAC-Bayes Inequality: Probabilistic and Information-Theoretic Foundations. The derivation of the PAC-Bayes inequality hinges on a combination of probabilistic tools and information-theoretic arguments. A central step involves applying a change of measure from $P$ to $Q$, leveraging the identity:

$$\mathbb{E}_{h \sim Q}[f(h)] = \mathbb{E}_{h \sim P} \left[ f(h) \frac{Q(h)}{P(h)} \right]. \qquad (4.7)$$
<figure><img src="image_38.png" /><figcaption>Figure 4.2: PAC-Bayes Bound vs. Sample Size (KL(Q||P) = 0.5)</figcaption></figure>

This allows the incorporation of $Q$ into bounds that originally apply to fixed $h$. By analyzing the moment-generating function of deviations between $\hat{R}(h, S)$ and $R(h)$, and applying Hoeffding's inequality to the empirical loss, we arrive at the following bound for any $Q$ and $P$, with probability at least $1 - \delta$:

$$R(Q) \le \hat{R}(Q, S) + \sqrt{\frac{\mathrm{KL}(Q||P) + \log \frac{1}{\delta}}{2m}}. \quad (4.8)$$

The generalization bound is therefore given by:

$$\mathcal{L}(f) - \mathcal{L}_{\mathrm{emp}}(f) \le \sqrt{\frac{\mathcal{KL}(Q||P) + \log(1/\delta)}{2N}}, \quad (4.9)$$

where $\mathcal{KL}(Q||P)$ quantifies the divergence between the posterior $Q$ and prior $P$. This bound is remarkable because it explicitly ties the true risk $R(Q)$ to the empirical risk $\hat{R}(Q, S)$, the KL divergence, and the sample size $m$. The PAC-Bayes bound encapsulates three competing forces: the empirical risk $\hat{R}(Q, S)$, the complexity penalty $\mathrm{KL}(Q||P)$, and the confidence term $\sqrt{\frac{\log \frac{1}{\delta}}{2m}}$. This interplay reflects a fundamental trade-off in learning:

1. **Empirical Risk:** $\hat{R}(Q, S)$ captures how well the posterior $Q$ fits the training data.

2. **Complexity:** The KL divergence ensures that $Q$ remains close to $P$, discouraging overfitting and promoting generalization.

3. **Confidence:** The term $\sqrt{\frac{\log \frac{1}{\delta}}{2m}}$ shrinks with increasing sample size, tightening the bound and enhancing reliability.
<figure><img src="image_39.png" /><figcaption>Figure 4.3: PAC-Bayes Bound vs. Sample Size (KL(Q||P) = 5.0)</figcaption></figure>

The KL term also introduces an inherent regularization effect, penalizing hypotheses that deviate significantly from prior knowledge. This aligns with Occam's Razor, favoring simpler explanations that are consistent with the data.

There are several extensions and Advanced Applications of Pac-Bayes Formalism. While the classical PAC-Bayes framework assumes i.i.d. data, recent advancements have generalized the theory to handle structured data, such as in time-series and graph-based learning. Furthermore, alternative divergence measures, like Rényi divergence or Wasserstein distance, have been explored to accommodate scenarios where KL divergence may be inappropriate. In practical settings, PAC-Bayes bounds have been instrumental in analyzing neural networks, Bayesian ensembles, and stochastic processes, offering theoretical guarantees even in high-dimensional, non-convex optimization landscapes.

### 4.1.2.1 Python Code to generate Figures 4.1, 4.2, and 4.3 illustrating PAC-Bayes Bound vs. Sample Size

The Python code below produces the Figures 4.1, 4.2, and 4.3 illustrating PAC-Bayes Bound vs. Sample Size.

```python
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimization import brentq

# --- KL divergence between Bernoulli(p) and Bernoulli(q) ---
def kl_bern(p, q):
```
```python
eps = 1e-12
p = np clipped(p, eps, 1-eps)
q = np clipped(q, eps, 1-eps)
return p * np log(p/q) + (1-p) * np log((1-p)/(1-q))

# --- PAC-Bayes bound inversion helper ---
def pac_bayes_bound(risk, kl_qp, n, delta=0.05):
    rhs = (kl_qp + np log(2* np sqrt(n)/delta)) / n

    # If RHS is too large, return trivial bound = 1
    if rhs >= np log(2):
        return 1.0

    # Solve for q >= risk s.t. KL(risk || q) = rhs
    def f(q):
        return kl_bern(risk, q) - rhs

    try:
        return brentq(f, risk+1e-9, 1-1e-9)
    except ValueError:
        return 1.0

# Parameters
delta = 0.05
emp_risks = np linspace(0.0, 0.5, 100)
kl_values = [0.5, 2.0, 5.0]
sample_sizes = [100, 500, 1000, 5000]

fig, axes = plt subplots(len(kl_values), 1, figsize=(8, 5*len(kl_values)),
    sharex = True)

for idx, kl in enumerate(kl_values):
    ax = axes[idx]
    for n in sample_sizes:
        bounds = [pac_bayes_bound(r, kl, n, delta) for r in emp_risks]
        ax plotzas risks, bounds, label=f"n={n}")
        ax set title(f"PAC-Bayes Bound (KL(Q||P)={kl})")
        ax set ylabel("True Risk Bound $L(Q)$")
        ax grid(True)
        ax legend()

axes[-1].set xlabel("Emprical Risk $\hat{L}(Q)$")
plt tight layout()
plt show()
```

### 4.1.3 KL divergence

The Kullback-Leibler (KL) divergence, also known as the relative entropy, is a fundamental concept in information theory, probability theory, and statistics. It quantifies the difference between two probability distributions, measuring the inefficiency of assuming that data is generated from one distribution when it is actually generated from another. Mathematically, for two discrete probability distributions $P(x)$ and $Q(x)$ defined over a common sample space $\mathcal{X}$, the KL divergence is
<figure><img src="image_40.png" /><figcaption>Figure 4.4: KL Divergence between Bernoulli(p) and Bernoulli(q)</figcaption></figure>

given by

$$D_{\text{KL}}(P \parallel Q) = \sum_{x \in \mathcal{X}} P(x) \log \frac{P(x)}{Q(x)} \qquad (4.10)$$

where the logarithm is conventionally taken to be the natural logarithm unless otherwise specified. For continuous probability distributions with probability density functions $p(x)$ and $q(x)$, the KL divergence is defined as

$$D_{\text{KL}}(P \parallel Q) = \int_{\mathcal{X}} p(x) \log \frac{p(x)}{q(x)} dx. \qquad (4.11)$$

The KL divergence can be understood as the expectation of the logarithmic difference between the two distributions, weighted by the true distribution $P(x)$:

$$D_{\text{KL}}(P \parallel Q) = \mathbb{E}_{P} \left[ \log \frac{P(x)}{Q(x)} \right]. \qquad (4.12)$$

Expanding this expectation in integral form,

$$D_{\text{KL}}(P \parallel Q) = \int_{\mathcal{X}} p(x) \log p(x) dx - \int_{\mathcal{X}} p(x) \log q(x) dx. \qquad (4.13)$$

The first integral represents the entropy of $P(x)$, denoted as

$$H(P) = - \int_{\mathcal{X}} p(x) \log p(x) dx, \qquad (4.14)$$
<figure><img src="image_41.png" /><figcaption>Figure 4.5: KL Divergence as function of q for fixed p</figcaption></figure>

while the second integral represents the cross-entropy between $P(x)$ and $Q(x)$, given by

$$H(P, Q) = - \int_{\mathcal{X}} p(x) \log q(x) dx. \qquad (4.15)$$

Thus, KL divergence can be rewritten as

$$D_{\mathrm{KL}}(P \parallel Q) = H(P, Q) - H(P), \qquad (4.16)$$

which represents the additional number of bits required to encode data from $P(x)$ when using a coding scheme optimized for $Q(x)$ instead of $P(x)$. Since the logarithm function is concave, an application of Jensen's inequality to the function $f(x) = -\log x$ yields the fundamental property

$$D_{\mathrm{KL}}(P \parallel Q) \geq 0, \qquad (4.17)$$

with equality if and only if $P(x) = Q(x)$ for all $x$, which follows from the strict convexity of the logarithm function. Unlike a metric, KL divergence is not symmetric, meaning that in general,

$$D_{\mathrm{KL}}(P \parallel Q) \neq D_{\mathrm{KL}}(Q \parallel P). \qquad (4.18)$$

For two multivariate normal distributions $P = N(\mu_1, \Sigma_1)$ and $Q = N(\mu_2, \Sigma_2)$, the KL divergence is given by

$$D_{\mathrm{KL}}(P \parallel Q) = \frac{1}{2} \left[ \mathrm{tr}(\Sigma_2^{-1} \Sigma_1) + (\mu_2 - \mu_1)^T \Sigma_2^{-1} (\mu_2 - \mu_1) - k + \log \frac{\det \Sigma_2}{\det \Sigma_1} \right], \qquad (4.19)$$
<figure><img src="image_42.png" /><figcaption>Figure 4.6: Gaussian KL Divergence (equal variances)</figcaption></figure>

where $k$ is the dimension of the Gaussian distribution. This expression reveals that KL divergence accounts for both the difference in mean vectors and the difference in covariance structures. The KL divergence has a fundamental connection to maximum likelihood estimation (MLE). Given a true data distribution $P(x)$ and a parametric model $Q(x | \theta)$, the parameter $\theta$ that minimizes the KL divergence

$$D_{\text{KL}}(P \parallel Q_{\theta}) = \mathbb{E}_{P}[\log P(x) - \log Q(x | \theta)] \quad (4.20)$$

is equivalent to the maximum likelihood estimate, since minimizing KL divergence is equivalent to maximizing the log-likelihood function

$$\sum_{i=1}^{N} \log Q(x_i | \theta). \quad (4.21)$$

This follows from the fact that the term $\mathbb{E}_{P}[\log P(x)]$ does not depend on $\theta$ and is therefore irrelevant for optimization. In Bayesian inference, KL divergence plays a central role in variational inference. Given an intractable posterior $P(z | x)$, one approximates it with a tractable distribution $Q(z)$, where the objective is to minimize

$$D_{\text{KL}}(Q(z) \parallel P(z | x)). \quad (4.22)$$

This minimization leads to the Evidence Lower Bound (ELBO),

$$\log P(x) \geq \mathbb{E}_{Q}[\log P(x, z)] - \mathbb{E}_{Q}[\log Q(z)], \quad (4.23)$$
<figure><img src="image_43.png" /><figcaption>Figure 4.7: Gaussian KL Divergence vs $\sigma_p$ ($\mu_q = \mu_p = 0$)</figcaption></figure>

which is central in modern probabilistic modeling. Moreover, KL divergence is related to mutual information. Given two random variables $X$ and $Y$ with a joint distribution $P(X, Y)$ and marginals $P(X)$ and $P(Y)$, the mutual information is given by

$$I(X; Y) = D_{\text{KL}}(P(X, Y) \parallel P(X)P(Y)), \qquad (4.24)$$

which quantifies the amount of information shared between $X$ and $Y$. In information geometry, the Fisher information matrix is derived from the second-order expansion of KL divergence,

$$g_{ij} = \mathbb{E}_P \left[ \frac{\partial \log P(x | \theta)}{\partial \theta^i} \frac{\partial \log P(x | \theta)}{\partial \theta^j} \right]. \qquad (4.25)$$

KL divergence is also connected to the Pinsker inequality,

$$\|P - Q\|_{\text{TV}} \le \sqrt{\frac{1}{2} D_{\text{KL}}(P \parallel Q)}. \qquad (4.26)$$

This establishes KL divergence as a fundamental measure of probability distance in statistical inference and machine learning.

The Kraft–McMillan inequality establishes a necessary and sufficient condition for the existence of a uniquely decodable code with a given set of codeword lengths. For a code over a D-ary alphabet to be uniquely decodable, the codeword lengths $l_1, l_2, ..., l_N$ must satisfy the inequality

$$\sum_{i=1}^{N} D^{-l_i} \le 1 \qquad (4.27)$$
This inequality is foundational because it directly links the lengths of codewords to a probability distribution. Specifically, if we define $m_i = D^{-l_i}$, then the Kraft sum condition

$$ \sum_{i} m_i \le 1 \qquad (4.28) $$

implies that $m_i$ can be normalized to form a probability distribution. This provides a formal bridge between the theory of codes and the theory of probability distributions, which is the central domain of information theory measures like the Kullback-Leibler (KL) Divergence.

The KL Divergence, or relative entropy, between two probability distributions $P$ and $Q$ over the same discrete set $\mathcal{X}$ is defined as

$$ D_{KL}(P \parallel Q) = \sum_{x \in \mathcal{X}} P(x) \log \frac{P(x)}{Q(x)} \qquad (4.29) $$

Its importance in data compression is given by the source coding theorem, which states that the expected codeword length $L(C)$ of an optimal code for a source with distribution $P$ is bounded by

$$ H(P) \le L(C) < H(P) + 1 \qquad (4.30) $$

where $H(P)$ is the Shannon entropy. When the code is designed for an incorrect distribution $Q$, the expected codeword length becomes

$$ L(C, Q) = \sum_{x} P(x)l_{x} \qquad (4.31) $$

The Kraft-McMillan theorem guarantees that for any uniquely decodable code with lengths $l_x$, there exists a probability distribution

$$ R(x) = D^{-l_x}/c \qquad (4.32) $$

where $c \le 1$ is the Kraft sum. This allows the lengths to be expressed as

$$ l_x = -\log_D R(x) + \log_D c \qquad (4.33) $$

and since $\log_D c \le 0$, we have

$$ l_x \ge -\log_D R(x) \qquad (4.34) $$

The expected length when using a code optimal for $Q$ on data from $P$ is

$$ L(C, Q) = \sum_{x} P(x)l_{x} \qquad (4.35) $$

Substituting the fundamental relation $l_x = -\log Q(x)$ for a code that meets the Kraft inequality with equality for distribution $Q$ yields

$$ L(C, Q) = - \sum_{x} P(x) \log Q(x) \qquad (4.36) $$

The difference between this suboptimal expected length and the entropy lower bound is

$$ L(C, Q) - H(P) = - \sum_{x} P(x) \log Q(x) + \sum_{x} P(x) \log P(x) \qquad (4.37) $$

This simplifies directly to $D_{KL}(P \parallel Q)$, as

$$ \sum_{x} P(x) \log \frac{P(x)}{Q(x)} = D_{KL}(P \parallel Q) \qquad (4.38) $$
Therefore, the KL Divergence quantifies the inefficiency, in bits or nats, incurred by using a code designed for a distribution $Q$ when the true data-generating distribution is $P$. The mathematical rigor of this relationship is entirely contingent upon the Kraft-McMillan theorem, which provides the necessary condition

$$ \sum_x D^{-l_x} \le 1 \qquad (4.39) $$

that legitimizes the interpretation of code lengths $l_x$ as representing a probability distribution, thereby enabling the derivation

$$ L(C, Q) = H(P) + D_{KL}(P \parallel Q) \qquad (4.40) $$

for a code that is optimal for $Q$.

### 4.1.3.1 Python Code to generate Figures 4.4, 4.5 illustrating KL Divergence in Bernoulli Random Variables and Figures 4.6, 4.7 illustrating KL Divergence in Gaussian Random Variables

The Python code below produces the Figures 4.4, 4.5 illustrating KL Divergence in Bernoulli Random Variables and Figures 4.6, 4.7 illustrating KL Divergence in Gaussian Random Variables.

```python
import numpy as np
import matplotlib.pyplot as plt

# --- Bernoulli KL divergence ---
def kl_bern(p, q):
    eps = 1e-12
    p = np clipped(p, eps, 1-eps)
    q = np clipped(q, eps, 1-eps)
    return p* np log(p/q) + (1-p) * np log((1-p)/(1-q))

# --- Gaussian KL divergence ---
def kl_gauss(mu_q, sigma_q, mu_p, sigma_p):
    return np logsigma_p/sigma_q) + (sigma_q**2 + (mu_q - mu_p)**2)/(2*sigma_p**2) - 0.5

# --- Bernoulli KL plots ---
ps = np linspace(0.001, 0.999, 100)
qs = np linspace(0.001, 0.999, 100)
P, Q = np meshgrid(ps, qs)
KL = kl_bern(P, Q)

fig, ax = plt figure(figsize=(7,6))
ax contourf(P, Q, KL, levels=50, cmap="viridis")
ax colorbar label="KL(Bern(p) || Bern(q))")
ax xlabel("p")
ax ylabel("q")
ax title("KL Divergence between Bernoulli(p) and Bernoulli(q)")
ax show()

fig, ax = plt figure(figsize=(8,6))
for p in [0.1, 0.3, 0.5, 0.7, 0.9]:
    ax plot(qs, kl_bern(p, qs), label=f"p={p}")
ax xlabel("q")
ax ylabel("KL(Bern(p) || Bern(q))")
```
```python
34 plt.title("KL Divergence as function of q for fixed p")
35 plt.legend()
36 plt.grid(True)
37 plt.show()
38
39 # --- Gaussian KL plots ---
40 mu_qq = np.linspace(-2, 2, 100)
41 mu_pq = np.linspace(-2, 2, 100)
42 MUQ, MUP = np.meshgrid(mu_qq, mu_pq)
43
44 # Fix sigma_q, sigma_p
45 sigma_q = 1.0
46 sigma_p = 1.0
47 KL_gauss = kl_gauss(MUQ, sigma_q, MUP, sigma_p)
48
49 plt.figure(figsize=(7,6))
50 plt.contourf(MUQ, MUP, KL_gauss, levels=50, cmap="plasma")
51 plt.colorbar label="KL(N($\mu_q$,1) || N($\mu_p$,1))")
52 plt.xlabel("$\mu_q$")
53 plt.ylabel("$\mu_p$")
54 plt.title("Gaussian KL Divergence (equal variances)")
55 plt.show()
56
57 # Vary sigma values
58 sigmas = [0.5, 1.0, 2.0]
59 plt.figure(figsize=(8,6))
60 mu_q = 0.0
61 mu_p = 0.0
62 for sigma_q in sigmas:
63     vals = [kl_gauss(mu_q, sigma_q, mu_p, s) for s in np.linspace(0.1,
64         3.0, 200)]
64     plt.plot(np.linspace(0.1, 3.0, 200), vals, label=f"$\sigma_q$={
65         sigma_q}")
66     plt.xlabel("$\sigma_p$")
67     plt.ylabel("KL$(N(\mu_q,\sigma_q^2) || N(\mu_p,\sigma_p^2))")
68     plt.title("Gaussian KL Divergence vs $\sigma_p$ ($\mu_q=\mu_p=0$)")
69     plt.legend()
70     plt.grid(True)
71     plt.show()
```

### 4.1.4 Rényi divergence

The Rényi divergence is a parametric family of divergences that generalizes the Kullback-Leibler (KL) divergence and provides a broader measure of dissimilarity between probability distributions.

Given two probability distributions $P$ and $Q$ defined on a measurable space $(\Omega, \mathcal{F})$, the Rényi divergence of order $\alpha$ (for $\alpha > 0$ and $\alpha \neq 1$) is defined as

$$D_{\alpha}(P \| Q) = \frac{1}{\alpha - 1} \ln \int_{\Omega} p(x)^{\alpha} q(x)^{1-\alpha} d\mu(x), \quad (4.41)$$

where $p(x)$ and $q(x)$ denote the Radon-Nikodym derivatives (i.e., probability density functions) of $P$ and $Q$ with respect to a common reference measure $\mu$. That is,

$$p(x) = \frac{dP}{d\mu}, \quad q(x) = \frac{dQ}{d\mu}. \qquad (4.42)$$
<figure><img src="image_44.png" /><figcaption>Figure 4.8: Rényi Divergence between Gaussians vs. α</figcaption></figure>

For the case where $P$ and $Q$ are discrete distributions supported on a finite or countable set $\mathcal{X}$, the Rényi divergence takes the form

$$D_{\alpha}(P \| Q) = \frac{1}{\alpha - 1} \ln \sum_{x \in \mathcal{X}} p(x)^{\alpha} q(x)^{1-\alpha}. \qquad (4.43)$$

Equivalently, it can be expressed using expectation notation as

$$D_{\alpha}(P \| Q) = \frac{1}{\alpha - 1} \ln \mathbb{E}_{Q} \left[ \left( \frac{p(X)}{q(X)} \right)^{\alpha} \right]. \qquad (4.44)$$

This expression shows that Rényi divergence is a moment-based measure of dissimilarity between $P$ and $Q$, generalizing the standard KL divergence, which is recovered in the limit $\alpha \to 1$. To see this, we differentiate

$$D_{\alpha}(P \| Q) = \frac{1}{\alpha - 1} \ln \int_{\Omega} p(x)^{\alpha} q(x)^{1-\alpha} d\mu(x), \qquad (4.45)$$

with respect to $\alpha$ and apply L'Hôpital's rule. Explicitly, defining

$$F(\alpha) = \int_{\Omega} p(x)^{\alpha} q(x)^{1-\alpha} d\mu(x), \qquad (4.46)$$

we compute

$$\frac{d}{d\alpha} F(\alpha) = \int_{\Omega} p(x)^{\alpha} q(x)^{1-\alpha} \ln \frac{p(x)}{q(x)} d\mu(x). \qquad (4.47)$$
<figure><img src="image_45.png" /><figcaption>Figure 4.9: Rényi Divergence Heatmap ($\alpha = 2.0, \sigma_p = \sigma_q = 1$)</figcaption></figure>

Applying the logarithmic derivative, we obtain

$$ \lim_{\alpha \to 1} D_{\alpha}(P\|Q) = \sum_{x \in \mathcal{X}} p(x) \ln \frac{p(x)}{q(x)} = D_{\text{KL}}(P\|Q), \quad (4.48) $$

demonstrating the convergence to the KL divergence. The Rényi divergence satisfies the following fundamental mathematical properties. It is non-negative,

$$ D_{\alpha}(P\|Q) \ge 0, \quad \text{with equality if and only if } P = Q. \quad (4.49) $$

Additionally, Rényi divergence is a monotonic function of $\alpha$, satisfying

$$ \alpha_1 < \alpha_2 \Rightarrow D_{\alpha_1}(P\|Q) \le D_{\alpha_2}(P\|Q). \quad (4.50) $$

For special values of $\alpha$, Rényi divergence simplifies to various known divergences. When $\alpha = 0$, it becomes

$$ D_0(P\|Q) = -\ln \sum_{\{x:p(x)>0\}} q(x), \quad (4.51) $$

which represents the negative logarithm of the probability mass that $Q$ assigns to the support of $P$. When $\alpha = 2$, we obtain

$$ D_2(P\|Q) = \ln \sum_{x \in \mathcal{X}} \frac{p(x)^2}{q(x)}, \quad (4.52) $$
which is related to the chi-squared divergence. When $\alpha = \infty$, the divergence simplifies to

$$D_{\infty}(P\|Q) = \ln \sup_{x \in \mathcal{X}} \frac{p(x)}{q(x)}, \qquad (4.53)$$

which measures the worst-case discrepancy between the two distributions. One crucial property of the Rényi divergence is its invariance under transformations. Given a Markov kernel $T$ that maps $P$ and $Q$ to new distributions $TP$ and $TQ$, the data-processing inequality states

$$D_{\alpha}(TP\|TQ) \leq D_{\alpha}(P\|Q), \qquad (4.54)$$

ensuring that no transformation can increase the divergence between distributions. Rényi divergence also satisfies the joint convexity property: for a mixture of probability distributions $P_i$ and $Q_i$ with weights $\lambda_i$,

$$D_{\alpha}\left(\sum_{i} \lambda_{i} P_{i} \middle\| \sum_{i} \lambda_{i} Q_{i}\right) \leq \sum_{i} \lambda_{i} D_{\alpha}(P_{i} \| Q_{i}). \qquad (4.55)$$

This inequality demonstrates that Rényi divergence behaves well under probabilistic mixing. Another fundamental result is the Pinsker-type bound, which states that for $\alpha > 1$, the Rényi divergence is bounded below by a function of the total variation distance $d_{TV}(P, Q)$,

$$D_{\alpha}(P\|Q) \geq \frac{\alpha}{\alpha - 1} \ln \left( 1 + \frac{\alpha - 1}{\alpha} d_{\text{TV}}(P, Q)^2 \right). \qquad (4.56)$$

For small $d_{TV}(P, Q)$, this inequality shows that Rényi divergence behaves quadratically in the deviation between $P$ and $Q$. Furthermore, the Rényi divergence is closely related to the Chernoff information, which determines the optimal exponent for Bayesian hypothesis testing. The Chernoff information is defined as

$$C(P, Q) = - \min_{\alpha \in (0,1)} D_{\alpha}(P \| Q), \qquad (4.57)$$

which plays a fundamental role in large deviations theory and hypothesis testing.

In summary, the Rényi divergence is a highly versatile and mathematically rich measure of dissimilarity between probability distributions. It generalizes the KL divergence, interpolates between various known divergence measures, satisfies essential properties such as monotonicity and data-processing inequalities, and has deep connections to large deviation theory, statistical estimation, and hypothesis testing.

### 4.1.4.1 Python Code to generate Figures 4.8 and 4.9 illustrating Rényi divergence

The Python code below produces the Figures 4.8 and 4.9 illustrating Rényi divergence.

```python
import numpy as np
import matplotlib.pyplot as plt

# --- Renyi divergence for univariate Gaussians ---
def renyi_gauss(mu_p, sigma_p, mu_q, sigma_q, alpha):
    if alpha == 1:
        raise ValueError("Alpha=1 gives KL divergence, not Renyi.")
    sigma alpha2 = alpha * sigma_p**2 + (1 - alpha) * sigma_q**2
    term1 = (alpha / (2*(1 - alpha))) * np.logsigma_q**2 / sigma alpha2)
    term2 = (alpha * (mu_p - mu_q)**2) / (2*sigma alpha2)
    return term1 + term2

# Parameters
```
14 alphas = np.linspace(0.1, 5.0, 200)
15 mu_p, mu_q = 0.0, 1.0
16 sigma_p, sigma_q = 1.0, 1.0
17
18 # --- Plot Renyi divergence vs alpha ---
19 vals = [renyi_gauss(mu_p, sigma_p, mu_q, sigma_q, a) for a in alphas]
20 plt.figure(figsize=(8,6))
21 plt.plot(alphas, vals)
22 plt.tight_layout()
23 plt.title("Renyi Divergence between Gaussians vs. $\alpha")
24 plt.grid(True)
25 plt.show()
26
27
28 # --- Heatmap over mu_p, mu_q for fixed alpha ---
29 mu vals = np.linspace(-2, 2, 100)
30 MU_P, MU_Q = np.meshgrid(mu vals, mu vals)
31 alpha = 2.0
32
33 REN = np zeros_like(MU_P)
34 for i in range(len(mu vals)):
35     for j in range(len(mu vals)):
36         REN[i,j] = renyi_gauss(MU_P[i,j], sigma_p, MU_Q[i,j], sigma_q, alpha)
37
38 plt figure(figsize=(7,6))
39 plt.contourf(MU_P, MU_Q, REN, levels=50, cmap="magma")
40 plt colorbar label=fr"Renyi Divergence $D_{\{\alpha\}}")"
41 plt.tight_layout()
42
43 # --- Plot Renyi divergence vs alpha ---
44

4.1.5 Wasserstein distance

The Wasserstein distance, also referred to as the Earth Mover's Distance (EMD), is a fundamental metric in the field of optimal transport theory that rigorously quantifies the distance between two probability measures by determining the minimal cost required to transport one measure into another.

Consider two probability distributions $\mu$ and $\nu$, defined on a metric space $(\mathcal{X}, d)$. The Wasserstein distance of order $p \ge 1$ is defined by

$$W_p(\mu, \nu) = \left( \inf_{\gamma \in \Pi(\mu, \nu)} \int_{\mathcal{X} \times \mathcal{X}} d(x, y)^p d\gamma(x, y) \right)^{\frac{1}{p}} \quad (4.58)$$

where $\Pi(\mu, \nu)$ denotes the set of all probability measures $\gamma$ on the product space $\mathcal{X} \times \mathcal{X}$ satisfying

$$\int_{\mathcal{X}} d\gamma(x, y) = \mu(x), \quad \int_{\mathcal{X}} d\gamma(x, y) = \nu(y). \qquad (4.59)$$

This definition arises from the optimal transport problem, originally formulated by Monge in 1781 and later extended by Kantorovich. The Monge formulation seeks a transport map $T: \mathcal{X} \to \mathcal{X}$
<figure><img src="image_46.png" /><figcaption>Figure 4.10: Wasserstein Distance vs Mean Difference ($\sigma_p = \sigma_q = 1$)</figcaption></figure>

that pushes $\mu$ onto $\nu$, i.e., satisfies

$$T_{\#}\mu = \nu \qquad (4.60)$$

where $T_{\#}\mu$ denotes the pushforward measure defined by

$$(T_{\#}\mu)(A) = \mu(T^{-1}(A)), \quad \forall A \subseteq \mathcal{X}. \qquad (4.61)$$

The Monge formulation minimizes the total transport cost given by

$$\inf_{T:T_{\#}\mu=\nu} \int_{\mathcal{X}} d(x, T(x))^p d\mu(x). \qquad (4.62)$$

However, the existence of an optimal transport map is highly constrained and may fail in many cases. To overcome this, Kantorovich introduced a relaxation in which transport maps are replaced by probability couplings $\gamma$, leading to the variational problem

$$W_p^p(\mu, \nu) = \inf_{\gamma \in \Pi(\mu, \nu)} \int_{\mathcal{X} \times \mathcal{X}} d(x, y)^p d\gamma(x, y). \qquad (4.63)$$

This formulation allows for mass splitting, making the problem more tractable and ensuring the existence of minimizers under mild conditions. The dual formulation of the Wasserstein distance is given by

$$W_p^p(\mu, \nu) = \sup_{\varphi, \psi} \left\{ \int_{\mathcal{X}} \varphi(x) d\mu(x) + \int_{\mathcal{X}} \psi(y) d\nu(y) \right\} \qquad (4.64)$$

subject to the constraint

$$\varphi(x) + \psi(y) \le d(x, y)^p, \quad \forall x, y \in \mathcal{X}. \qquad (4.65)$$
<figure><img src="image_47.png" /><figcaption>Figure 4.11: Wasserstein Distance vs Variance ($\mu_p = \mu_q = 0, \sigma_q = 1$)</figcaption></figure>

For $p = 1$, the Kantorovich-Rubinstein theorem provides an alternative characterization

$$W_1(\mu, \nu) = \sup_{\|\nabla \varphi\| \le 1} \left\{ \int_X \varphi(x) d\mu(x) - \int_X \varphi(x) d\nu(x) \right\}. \qquad (4.66)$$

When $\mathcal{X} = \mathbb{R}$, the Wasserstein-1 distance has the closed-form expression

$$W_1(\mu, \nu) = \int_{-\infty}^{\infty} |F_{\mu}(x) - F_{\nu}(x)| dx \qquad (4.67)$$

where $F_\mu$ and $F_\nu$ are the cumulative distribution functions of $\mu$ and $\nu$, respectively. For $p = 2$, the squared Wasserstein-2 distance is given by

$$W_2^2(\mu, \nu) = \inf_{\gamma \in \Pi(\mu, \nu)} \int_{\mathcal{X} \times \mathcal{X}} d(x, y)^2 d\gamma(x, y). \qquad (4.68)$$

If $\mu$ and $\nu$ are Gaussian distributions $\mathcal{N}(m_1, \Sigma_1)$ and $\mathcal{N}(m_2, \Sigma_2)$, then

$$W_2^2(\mu, \nu) = \|m_1 - m_2\|^2 + \mathrm{Tr}(\Sigma_1 + \Sigma_2 - 2(\Sigma_1^{1/2} \Sigma_2 \Sigma_1^{1/2})^{1/2}). \qquad (4.69)$$

The Wasserstein distance induces a metric on the space of probability measures, providing a notion of distance that is weaker than total variation but stronger than weak convergence. Specifically, if $\mu_n \to \mu$ in Wasserstein distance, then for any function $f$ satisfying

$$|f(x)| \le C(1 + d(x, x_0)^p), \qquad (4.70)$$
<figure><img src="image_48.png" /><figcaption>Figure 4.12: Wasserstein Distance Heatmap ($\sigma_p = 1.0, \sigma_q = 1.5$)</figcaption></figure>

it holds that

$$ \int_X f(x)d\mu_n(x) \rightarrow \int_X f(x)d\mu(x). \qquad (4.71) $$

The Wasserstein metric endows the space of probability measures with a Riemannian-like structure, where geodesics are given by **McCann interpolation**, defined as

$$ \mu_t = ((1 - t) \mathrm{id} + tT)_\# \mu_0, \quad t \in [0, 1]. \qquad (4.72) $$

For discrete distributions

$$ \mu = \sum_{i=1}^{n} a_i \delta_{x_i}, \quad \nu = \sum_{j=1}^{m} b_j \delta_{y_j}, \qquad (4.73) $$

the Wasserstein distance reduces to a linear programming problem

$$ \min_{\gamma \in \Pi(\mu, \nu)} \sum_{i,j} \gamma_{ij} d(x_i, y_j)^p. \qquad (4.74) $$

Computationally, Wasserstein distances are expensive to evaluate, but regularized approximations such as **Sinkhorn's algorithm** solve the entropically regularized problem

$$ W_p^\lambda(\mu, \nu) = \inf_{\gamma \in \Pi(\mu, \nu)} \sum_{i,j} \gamma_{ij} d(x_i, y_j)^p + \lambda \sum_{i,j} \gamma_{ij} \log \gamma_{ij}. \qquad (4.75) $$

This approximation allows for efficient computation via iterative updates using the Sinkhorn-Knopp algorithm. The Wasserstein distance plays a crucial role in probability theory, statistics, machine
46 plt.title("Wasserstein Distance Heatmap ($\sigma_p$={sigma_p}, $\sigma_q$={sigma_q})")
47 plt.show()

## 4.2 Spectral Regularization

The concept of spectral regularization, which refers to the preferential learning of low-frequency modes by neural networks before high-frequency modes, emerges from a combination of Fourier analysis, optimization theory, and the inherent properties of deep neural networks. This phenomenon is tightly connected to the functional approximation capabilities of neural networks and can be rigorously understood through the lens of Fourier decomposition and the gradient descent optimization process.

### 4.2.1 Literature Review of Spectral Regularization

Jin et. al. (2025) [110] introduced a novel confusional spectral regularization technique to improve fairness in machine learning models. The study focuses on the spectral norm of the robust confusion matrix and proposes a method to control spectral properties, ensuring more robust and unbiased learning. It provides insights into how regularization can mitigate biases in classification tasks. Ye et. al. (2025) [111] applied spectral clustering with regularization to detect small clusters in complex networks. The work enhances spectral clustering techniques by integrating regularization methods, allowing improved performance in anomaly detection and community detection tasks. The approach significantly improves robustness in highly noisy data environments. Bhattacharjee and Bharadwaj (2025) [112] explored how spectral domain representations can benefit from autoencoder-based feature extraction combined with stochastic regularization techniques. The authors propose a Symmetric Autoencoder (SymAE) that enables better generalization of spectral features, particularly useful in high-dimensional data and deep learning applications. Wu et. al. (2025) [113] applied spectral regularization to geophysical data processing, specifically for high-resolution velocity spectrum analysis. The approach enhances the resolution of velocity estimation in seismic imaging by using hyperbolic Radon transform regularization, demonstrating how spectral regularization can benefit applications beyond traditional ML. Ortega et. al. (2025) [114] applied Tikhonov regularization to atmospheric spectral analysis, optimizing gas retrieval strategies in high-resolution spectroscopic observations. The work significantly improves methane (CH4) and nitrous oxide (N2O) detection accuracy by reducing noise in spectral measurements, showcasing the impact of spectral regularization in remote sensing and environmental monitoring. Kazmi et. al. (2025) [115] proposed a spectral regularization-based federated learning model to improve robustness in cybersecurity threat detection. The model addresses the issue of non-IID data in SDN (Software Defined Networks) by utilizing spectral norm-based regularization within deep learning architectures. Zhao et. al. (2025) [116] introduced a regularized deep spectral clustering method, which enhances feature selection and clustering robustness. The authors utilize projected adaptive feature selection combined with spectral graph regularization, improving clustering accuracy and interpretability in high-dimensional datasets. Saranya and Menaka (2025) [117] integrated spectral regularization with quantum-based machine learning to analyze EEG signals for Autism Spectrum Disorder (ASD) detection. The proposed method improves spatial filtering and feature extraction using wavelet-based regularization, leading to more reliable EEG pattern recognition. Cox and Ghosh (2022) [72] develop a partial factorization algorithm that utilizes Mertens' function to factor certain natural numbers, introduce a modified variant of Mertens' function constrained to negative values (except at $M(1)$), and explore its relationship with the sum-of-divisors function, thereby revealing new arithmetic connections between multiplicative and additive number-theoretic structures. Dhalbisoi et. al. (2024) [118] developed a Regularized Zero-Forcing (RZF) method for
spectral efficiency optimization in beyond 5G networks. The authors demonstrate that spectral regularization techniques can significantly improve signal-to-noise ratios in wireless communication systems, optimizing data transmission in massive MIMO architectures. Wei et. al. (2025) [119] explored the use of spectral regularization in medical imaging, particularly in 3D near-infrared spectral tomography. The proposed model integrates regularized convolutional neural networks (CNNs) to improve tissue imaging resolution and accuracy, demonstrating an application of spectral regularization in biomedical engineering.

### 4.2.2 Analysis of Spectral Regularization

Let us define a target function $f(\mathbf{x})$, where $\mathbf{x} \in \mathbb{R}^d$, and its Fourier transform $\hat{f}(\xi)$ as

$$ \hat{f}(\xi) = \int_{\mathbb{R}^d} f(\mathbf{x}) e^{-i2\pi\xi \cdot \mathbf{x}} d\mathbf{x} \qquad (4.76) $$

This transform breaks down $f(\mathbf{x})$ into frequency components indexed by $\xi$. In the context of deep learning, we seek to approximate $f(\mathbf{x})$ with a neural network output $f_{\text{NN}}(\mathbf{x}; \boldsymbol{\theta})$, where $\boldsymbol{\theta}$ represents the set of trainable parameters. The loss function to be minimized is typically the mean squared error:

$$ \mathcal{L}(\boldsymbol{\theta}) = \int_{\mathbb{R}^d} |f(\mathbf{x}) - f_{\text{NN}}(\mathbf{x}; \boldsymbol{\theta})|^2 d\mathbf{x} \qquad (4.77) $$

We can equivalently express this loss in the Fourier domain, leveraging Parseval's theorem:

$$ \mathcal{L}(\boldsymbol{\theta}) = \int_{\mathbb{R}^d} \left| \hat{f}(\boldsymbol{\xi}) - \hat{f}_{\text{NN}}(\boldsymbol{\xi}; \boldsymbol{\theta}) \right|^2 d\boldsymbol{\xi} \qquad (4.78) $$

To solve for $\boldsymbol{\theta}$, we employ gradient descent:

$$ \boldsymbol{\theta}^{(t+1)} = \boldsymbol{\theta}^{(t)} - \eta \nabla_{\boldsymbol{\theta}} \mathcal{L}(\boldsymbol{\theta}) \qquad (4.79) $$

where $\eta$ is the learning rate. The gradient of the loss function with respect to $\boldsymbol{\theta}$ is

$$ \nabla_{\boldsymbol{\theta}} \mathcal{L}(\boldsymbol{\theta}) = 2 \int_{\mathbb{R}^d} (\hat{f}_{\text{NN}}(\boldsymbol{\xi}; \boldsymbol{\theta}) - \hat{f}(\boldsymbol{\xi})) \nabla_{\boldsymbol{\theta}} \hat{f}_{\text{NN}}(\boldsymbol{\xi}; \boldsymbol{\theta}) d\boldsymbol{\xi} \qquad (4.80) $$

At the core of this gradient descent process lies the behavior of the gradient $\nabla_{\boldsymbol{\theta}} \hat{f}_{\text{NN}}(\boldsymbol{\xi}; \boldsymbol{\theta})$ with respect to the frequency components $\boldsymbol{\xi}$. For neural networks, particularly those with ReLU activations, the gradients of the output with respect to the parameters are expected to decay for high-frequency components. This can be approximated as

$$ \mathcal{R}(\xi) \sim \frac{1}{1 + \|\xi\|^2} \qquad (4.81) $$

which implies that the neural network is inherently more sensitive to low-frequency components of the target function during early iterations of training. This spectral decay is a direct consequence of the structure of the network's activations, which are more sensitive to low-frequency features due to their smoother, lower-order terms. To understand the role of the neural tangent kernel (NTK), which governs the linearized dynamics of the neural network, we define the NTK as

$$ \Theta(\mathbf{x}, \mathbf{x}'; \boldsymbol{\theta}) = \sum_{i=1}^{P} \frac{\partial f_{\text{NN}}(\mathbf{x}; \boldsymbol{\theta})}{\partial \theta_i} \frac{\partial f_{\text{NN}}(\mathbf{x}'; \boldsymbol{\theta})}{\partial \theta_i} \qquad (4.82) $$

The NTK essentially describes how the output of the network changes with respect to its parameters. The evolution of the network's output during training can be approximated by the solution to a linear system governed by the NTK. The output of the network at time $t$ is given by

$$ f_{\text{NN}}(\mathbf{x}; t) = \sum_{k} c_k (1 - e^{-\eta \lambda_k t}) \phi_k(\mathbf{x}) \qquad (4.83) $$
<figure><img src="image_49.png" /><figcaption>Figure 4.13: Effect of Spectral Regularization on Singular Values</figcaption></figure>

<figure><img src="image_50.png" /><figcaption>Figure 4.14: Original Weight Matrix & After Spectral Regularization</figcaption></figure>

where {$\lambda_k$} are the eigenvalues of $\Theta$, and {$\phi_k(x)$} are the corresponding eigenfunctions. The eigenvalues $\lambda_k$ determine the speed of convergence for each frequency mode, with low-frequency modes (large $\lambda_k$) converging more quickly than high-frequency ones (small $\lambda_k$):

$$1 - e^{-\eta\lambda_k t} \to 1 \quad \text{for large } \lambda_k \quad \text{and} \quad 1 - e^{-\eta\lambda_k t} \to 0 \quad \text{for small } \lambda_k \qquad (4.84)$$