(Continued from previous page)

<table><thead><tr><td>Authors (Year)</td><td>Contribution</td></tr></thead><tbody><tr><td>Lalor & Swishchuk (2025)</td><td>Extend SAC to non-Markovian market-making, demonstrating its effectiveness in handling long-term dependencies and stochastic pricing models.</td></tr><tr><td>Zhang et al. (2025)</td><td>Propose a diffusion-based SAC framework for multi-UAV networks in the Metaverse, optimizing cooperative task allocation and resource distribution.</td></tr><tr><td>Zhao et al. (2025)</td><td>Utilize SAC for energy management in hybrid storage systems for urban rail transit, optimizing power distribution in traction power supply systems.</td></tr><tr><td>Tresca et al. (2025)</td><td>Apply SAC to develop adaptive energy management strategies for hybrid electric vehicles, enhancing fuel efficiency and battery longevity through dynamic energy consumption adjustments.</td></tr></tbody></table>

### 16.4.6.7.3 Mathematical Analysis of Soft Actor-Critic (SAC) algorithm

The **Soft Actor-Critic (SAC)** algorithm is a model-free, off-policy deep reinforcement learning method that optimizes policies in **continuous action spaces**. It is formulated within the **maximum entropy reinforcement learning framework**, which aims to **maximize both the expected return and the entropy** of the policy. Given a Markov Decision Process (MDP) defined by the tuple $\mathcal{M} = (\mathcal{S}, \mathcal{A}, p, r, \gamma)$, where $\mathcal{S}$ is the state space, $\mathcal{A}$ is the action space, $p(s' | s, a)$ is the transition probability, $r(s, a)$ is the reward function, and $\gamma \in (0, 1]$ is the discount factor, SAC introduces a policy optimization objective that includes an entropy term. This results in a policy that is **stochastic** and encourages **exploration** while still optimizing for high rewards. The objective function for SAC is defined as

$$J(\pi) = \sum_{t=0}^{\infty} \mathbb{E}_{(s_t, a_t) \sim \rho_{\pi}} [r(s_t, a_t) + \alpha \mathcal{H}(\pi(\cdot | s_t))] \quad (16.340)$$

where $\rho_{\pi}(s_t, a_t)$ represents the state-action distribution under policy $\pi$, $\mathcal{H}(\pi(\cdot | s_t))$ is the entropy of the policy at state $s_t$, and $\alpha$ is a temperature parameter that controls the trade-off between **exploration and exploitation**. The soft Q-function in SAC satisfies the Bellman equation

$$Q^{\pi}(s, a) = r(s, a) + \gamma \mathbb{E}_{s' \sim p(s'|s, a), a' \sim \pi} [Q^{\pi}(s', a') - \alpha \log \pi(a' | s')] \quad (16.341)$$

which differs from the standard Q-function by incorporating the entropy term $-\alpha \log \pi(a' | s')$, making the policy **more stochastic** and **less greedy**. The policy is optimized by minimizing the KL divergence between $\pi(a | s)$ and an exponential Boltzmann distribution induced by $Q^{\pi}(s, a)$:

$$\pi^*(a | s) = \arg \min_{\pi} D_{\text{KL}} \left( \pi(\cdot | s) \middle\| \frac{\exp(Q(s, a)/\alpha)}{Z(s)} \right) \quad (16.342)$$

where $Z(s)$ is the partition function ensuring normalization. This results in a policy of the form

$$\pi(a | s) = \frac{\exp(Q(s, a)/\alpha)}{\int_{\mathcal{A}} \exp(Q(s, a')/\alpha) da'} \quad (16.343)$$

which highlights the **softmax-like behavior** of the policy. To **stabilize learning**, SAC employs **two Q-networks**, $Q_{\theta_1}(s, a)$ and $Q_{\theta_2}(s, a)$, and uses the minimum value in the Bellman update to reduce overestimation bias:

$$y = r(s, a) + \gamma \mathbb{E}_{s' \sim p, a' \sim \pi} [\min(Q_{\theta_1}(s', a'), Q_{\theta_2}(s', a')) - \alpha \log \pi(a' | s')] \quad (16.344)$$
<figure><img src="image_210.png" /><figcaption>Figure 16.42: Architecture of the Soft Actor-Critic (SAC) algorithm: The actor (green) generates actions from states, guided by entropy regularization (red) to encourage exploration. The critics (blue) estimate Q-values, updated using targets (yellow) for stability. A replay buffer (violet) stores past transitions from the environment (gray). Edges denote information flow, with labels indicating the quantities passed between components</figcaption></figure>

The loss function for training the Q-networks is given by

$$J_Q(\theta_i) = \mathbb{E}_{(s,a,r,s') \sim \mathcal{D}} \left[ (Q_{\theta_i}(s, a) - y)^2 \right] \quad (16.345)$$

where $\mathcal{D}$ is the **replay buffer**. The policy $\pi_{\phi}(a | s)$ is parameterized using a **neural network**, and its objective function is derived from the reparameterization trick. If $\pi_{\phi}$ is a Gaussian policy where

$$a = \tanh(\mu_{\phi}(s) + \sigma_{\phi}(s) \cdot \epsilon), \quad \epsilon \sim \mathcal{N}(0, I) \quad (16.346)$$

the policy loss is

$$J_{\pi}(\phi) = \mathbb{E}_{s \sim \mathcal{D}, \epsilon \sim \mathcal{N}(0, I)} \left[ \alpha \log \pi_{\phi}(a | s) - \min(Q_{\theta_1}(s, a), Q_{\theta_2}(s, a)) \right] \quad (16.347)$$

where the first term encourages **higher entropy**, and the second term ensures **policy improvement**. The temperature parameter $\alpha$ is **automatically adjusted** by minimizing

$$J_{\alpha}(\alpha) = \mathbb{E}_{s \sim \mathcal{D}, a \sim \pi_{\phi}} \left[ -\alpha \left( \log \pi_{\phi}(a | s) + \bar{H} \right) \right] \quad (16.348)$$

where $\bar{H}$ is a target entropy value. The gradients for $\alpha$ update are given by

$$\alpha \leftarrow \alpha - \lambda \nabla_{\alpha} J_{\alpha}(\alpha) \quad (16.349)$$
where $\lambda$ is the learning rate. The full SAC algorithm follows an **off-policy training procedure** where experience tuples $(s, a, r, s')$ are sampled from the replay buffer, and gradients are computed using stochastic updates. The key components include the **Q-function updates**, **policy updates**, and **entropy adjustments**, ensuring **stable convergence** and **exploratory behavior** throughout training.

### 16.4.6.7.4 Python Code to generate Figure 16.42

Soft Actor-Critic (SAC) is a reinforcement learning algorithm that can be visualized as a computational graph of modules rather than as a population of networks. A good way to plot SAC is to show the information flow between:

1. Actor (Policy Network) → outputs actions.

2. Critics (Q-networks) → estimate Q-values.

3. Target networks → stabilize learning.

4. Entropy term → adds exploration via temperature parameter $\alpha$

5. Replay buffer → stores experiences.

The Python code below produces the Figure 16.42 illustrating the Soft Actor-Critic (SAC).

```python
import matplotlib.pyplot as plt
import networkx as nx

def build_sac_graph():
    G = nx.DiGraph()

    # Nodes
    G.add_node("Replay Buffer", type="memory")
    G.add_node("Actor (Policy $\pi$)$", type="actor")
    G.add_node("Critic Q1$\Phi$", type="critic")
    G.add_node("Critic Q2$\Phi$", type="critic")
    G.add_node("Target Q1$\Phi$", type="target")
    G.add_node("Target Q2$\Phi$", type="target")
    G.add_node("Environment", type="env")
    G.add_node("Entropy $\alpha$", type="entropy")

    # Edges
    G.add_edge("Environment", "Replay Buffer", label=(s,a,r,s))
    G.add_edge("Replay Buffer", "Actor (Policy $\pi$)$", label="s")
    G.add_edge("Actor (Policy $\pi$)$", "Environment", label="a")

    G.add_edge("Replay Buffer", "Critic Q1$\Phi$", label=(s,a))
    G.add_edge("Replay Buffer", "Critic Q2$\Phi$", label=(s,a))

    G.add_edge("Actor (Policy $\pi$)$", "Critic Q1$\Phi$", label="a ~ $\pi$")
    G.add_edge("Actor (Policy $\pi$)$", "Critic Q2$\Phi$", label="a ~ $\pi$")

    G.add_edge("Target Q1$\Phi$", "Critic Q1$\Phi$", label="target")
    G.add_edge("Target Q2$\Phi$", "Critic Q2$\Phi$", label="target")

    G.add_edge("Entropy $\alpha$", "Actor (Policy $\pi$)$", label="exploration")
```
return G

def plot_sac save_as = None

G = build_sac_graph()

pos = {
    "Environment": (0,2),
    "Replay Buffer": (0,1),
    "Actor (Policy πθ)": (2,2),
    "Entropy α": (3,2.8),
    "Critic Q1Φ": (2,1),
    "Critic Q2Φ": (2,0.5),
    "Target Q1Φ": (4,1),
    "Target Q2Φ": (4,0.5),
}

# Colors

color_map = []

for n, d in G.nodes(data=True):
    if d["type"] == "actor":
        color_map.append("lightgreen")
    elif d["type"] == "critic":
        color_map.append("lightblue")
    elif d["type"] == "target":
        color_map.append("lightyellow")
    elif d["type"] == "entropy":
        color_map.append("salmon")
    elif d["type"] == "memory":
        color_map.append("violet")
    elif d["type"] == "env":
        color_map.append("lightgrey")
    else:
        color_map.append("white")

# Draw nodes

nx.draw_networkx_nodes(G, pos, node_color= color_map, edgecolors="k",
                       node_size=2500)

nx.draw_networkx labels (G, pos, font_size =9)

# Draw edges

ax = plt.gca()

for u, v, d in G. edges(data = True):
    ax. annotate("",
               xy = pos [v], xycoordinates = "data",
               xytext = pos [u], textcoordinates = "data",
               arrowprops = dict(arrowstyle = "->", color = "black", lw
                                 =1.5))

# Edge labels

if "label" in d:
    mid_x = (pos[u][0] + pos[v][0]) / 2
    mid_y = (pos[u][1] + pos[v][1]) / 2
    plt.text(mid_x, mid_y + 0.1, d["label"], fontsize = 7, ha = "center")

plt. title("Soft Actor-Critic (SAC) Architecture", fontsize =12)

plt. axis("off")
```
85
86
87
88
89
90
91

if save_as:
    plt.savefig save_as, dpi=300, bbox_inches="tight")
    plt.show()

if __name__ == "__main__":
    plot_sac save_as="sac_architecture.png")
```

## 16.5 Neuroevolution

Neuroevolution is a computational framework in artificial intelligence (AI) that utilizes evolutionary algorithms to optimize the architecture and weights of artificial neural networks (ANNs). Unlike conventional deep learning, which relies on gradient-based optimization techniques such as stochastic gradient descent (SGD), neuroevolution leverages principles of evolutionary computation, including mutation, crossover, and selection, to iteratively refine the structure and parameters of neural networks. Given a population of candidate neural networks $N = \{N_1, N_2, \dots, N_m\}$, each network is evaluated by a fitness function $F(N)$, which measures its performance on a given task. The fundamental process can be mathematically described as an iterative search over the space of possible networks, where the objective is to maximize the fitness function:

$$N^* = \arg \max_{N \in \mathcal{N}} F(N). \qquad (16.350)$$

Each neural network in the population is parameterized by a set of weights $W$ and biases $B$, where the function represented by the network is given by

$$y = \sigma(Wx + B), \qquad (16.351)$$

where $x$ is the input vector, $y$ is the output, and $\sigma(\cdot)$ is the activation function, typically chosen as a nonlinear function such as the sigmoid, tanh, or ReLU. The evolutionary algorithm modifies the parameters $W$ and $B$ over successive generations through genetic operators, which can be expressed as follows. Given a set of parent solutions $\{(W_i, B_i)\}_{i=1}^m$, offspring solutions are generated via mutation and crossover. Mutation introduces perturbations in the network weights:

$$W' = W + \epsilon, \quad B' = B + \delta, \qquad (16.352)$$

where $\epsilon \sim \mathcal{N}(0, \sigma^2)$ and $\delta \sim \mathcal{N}(0, \sigma^2)$ are Gaussian perturbations. The crossover operation combines weight matrices from two parent networks $(W_1, B_1)$ and $(W_2, B_2)$, producing an offspring $(W', B')$:

$$W' = \alpha W_1 + (1 - \alpha) W_2, \quad B' = \alpha B_1 + (1 - \alpha) B_2, \qquad (16.353)$$

where $\alpha \in [0, 1]$ is a crossover coefficient that determines the proportion of contribution from each parent. The fitness of each offspring network is then computed, and a selection mechanism is applied to determine which networks survive to the next generation. A common selection strategy is tournament selection, in which $k$ networks are randomly chosen, and the one with the highest fitness is selected:

$$N_{\text{next}} = \arg \max_{N \in \mathcal{N}_{\text{tournament}}} F(N). \qquad (16.354)$$

In addition to optimizing the weights, neuroevolution can be extended to evolve the architecture of the neural networks, including the number of layers, neurons per layer, and connectivity patterns. This results in a combinatorial search problem where each network is represented as a graph $G = (V, E)$, where $V$ is the set of neurons and $E$ is the set of synaptic connections. The evolution of architectures can be modeled using genetic encoding schemes such as direct encoding, where each network parameter is explicitly represented in a genome, or indirect encoding, where a compact
Neuroevolution: Parents Evolving into Offspring

<figure><img src="image_211.png" /><figcaption>Figure 16.43: Illustration of Neuroevolution: The top row shows neural network individuals from the parent generation, and the bottom row shows their offspring after evolutionary operations such as selection, crossover, and mutation. Nodes are color-coded by type (blue = input, green = hidden, red = output). Edge colors denote weight sign (blue = positive, red = negative), and thickness encodes weight magnitude</figcaption></figure>

representation (e.g., a developmental rule) is used to generate the architecture. If the architecture is encoded as a vector $\theta$, then the objective function is

$$ \theta^* = \arg \max_{\theta} F(G(\theta)). \qquad (16.355) $$

A powerful variant of neuroevolution is the use of novelty search, in which selection is based not on performance but on behavioral diversity. The novelty of a network is defined as the average distance between its behavior and those of its nearest neighbors in the population:

$$ \mathcal{N}(N) = \frac{1}{k} \sum_{i=1}^{k} d(N, N_i), \qquad (16.356) $$

where $d(N, N_i)$ is a distance metric (e.g., Euclidean or cosine distance) between network behaviors. Modern neuroevolution approaches incorporate gradient-based methods to accelerate convergence. One such method is Evolution Strategies (ES), which approximates the gradient of the expected fitness with respect to the network parameters via a perturbation-based estimation:

$$ \nabla_{\theta} \mathbb{E}[F(\theta)] \approx \frac{1}{\sigma} \sum_{i=1}^{n} F(\theta + \sigma \epsilon_i) \epsilon_i, \qquad (16.357) $$

where $\epsilon_i \sim \mathcal{N}(0, I)$ are random perturbations, and $\sigma$ is the step size. This update rule enables evolution to operate in high-dimensional parameter spaces more efficiently. Another advanced neuroevolution technique is the NeuroEvolution of Augmenting Topologies (NEAT), which evolves both network weights and topologies by introducing genetic operators such as mutation of connections and nodes. The key idea is to preserve structural innovations via speciation, where networks are clustered into species based on a similarity metric $S(N_1, N_2)$, typically computed as a function of the number of disjoint and excess connections:

$$ S(N_1, N_2) = c_1 E + c_2 D + c_3 W, \qquad (16.358) $$
where *E* and *D* are the number of excess and disjoint genes, respectively, and *W* is the average weight difference between matching genes, with *c*₁, *c*₂, *c*₃ as weighting coefficients. Through iterative application of these principles, neuroevolution generates increasingly effective neural networks for complex tasks such as reinforcement learning, robotics, and generative modeling. By searching the space of neural architectures and parameters simultaneously, it provides an alternative to traditional deep learning methods, enabling the discovery of novel architectures without the need for explicit human-designed features. The optimization process can be formally represented as

$$ \min_{\theta} \mathbb{E}_{\pi_{\theta}} \left[ \sum_{t=0}^{T} \gamma^t r_t \right], \qquad (16.359) $$

where $\pi_{\theta}$ is a policy parameterized by a neural network, $r_t$ is the reward at time step $t$, and $\gamma$ is a discount factor in reinforcement learning settings.

### 16.5.1 Python Code to generate Figure 16.43

Neuroevolution is the broad paradigm of evolving neural networks with evolutionary algorithms. A clear visualization usually has two parts:

1. Population of neural networks (at a given generation).

2. Evolutionary process: selection → crossover/mutation → next generation.

To keep it simple but illustrative, let's plot:

1. A population of small neural networks.

2. Arrows showing how parents produce offspring.

The Python code below produces the Figure 16.43 illustrating the Neuroevolution.

```python
import matplotlib.pyplot as plt
import networkx as nx
import random
```

```python
# ------------------------
# Build a simple neural network
# ------------------------
```

```python
def build_network(nInputs=2, nHidden=2, nOutputs=1):
    G = nx.DiGraph()
    # Input nodes
    for i in range(nInputs):
        G.add_node(f"I{i+1}", layer="input", pos=(i, 0))
    # Hidden nodes
    for i in range(nHidden):
        G.add_node(f"H{i+1}", layer="hidden", pos=(i, 1))
    # Output nodes
    for i in range(nOutputs):
        G.add_node(f"O{i+1}", layer="output", pos=(i, 2))
    # Random connections
    nodes = list(G.nodes())
    for _ in range(randomonde(3, 6)):
```
```python
u, v = random.sample nodes, 2)
if u != v:
    G.add_edge(u, v, weight = random .uniform(-1, 1))

return G
```
```python
draw_network(G, axes[1, i], title=f"Offspring {i+1}")

# Evolution arrows between generations
for i in range(3):
    axes[0, i]. annotate("",
                       xy=(0.5, -0.1), xycoordinates="axes fraction",
                       xytext=(0.5, -1.0), textcoordinates="axes
                                  fraction",
                       arrowprops=dict(arrowstyle="->", color="
                                                   black", lw=2))

plt.suptitle("Neuroevolution: Parents Evolving into Offspring",
             fontsize=14)
plt.tight_layout rect=[0, 0.03, 1, 0.95])

if save_as:
    plt.savefig save_as, dπ=300, bbox_inches="tight")
    plt.show()

# Example run
if __name__ == "__main__":
    plot_neuroevolution save_as="neuroevolution.png")
```

## 16.5.2 Neuro-Genetic Evolution

### 16.5.2.1 Literature Review of Neuro-Genetic Evolution

Neuro-genetic evolution, or neuroevolution, is a paradigm that integrates evolutionary algorithms with deep learning to optimize neural network architectures and training methodologies. The seminal work **NeuroEvolution of Augmenting Topologies (NEAT)** by Stanley and Miikkulainen (2002) [983] and Stanley et. al. (2005) [984] laid the foundation for evolving both the structure and weights of neural networks. NEAT employs a genetic algorithm that begins with minimal structures and incrementally complexifies them over generations while preserving novel architectures via speciation. Building on this, Gauci and Stanley (2007) [985] introduced **HyperNEAT** which is a compositional pattern-producing network (CPPN) to encode connectivity patterns, leading to the evolution of large-scale architectures with inherent symmetries and regularities. These methodologies significantly improved the scalability and efficiency of neuroevolution, inspiring subsequent research into evolving deep networks. **EANT (Evolutionary Acquisition of Neural Topologies)** by Kassahun and Sommer (2005) [987] took a similar approach, where networks start simple and gain complexity over time, thereby reducing the computational burden of evolving highly intricate architectures.

As deep learning gained prominence, researchers sought ways to automate the design of neural architectures. Miikkulainen et al. (2024) [991] introduced **CoDeepNEAT**, an extension of NEAT that optimizes deep learning architectures by evolving topology, hyperparameters, and components, yielding results comparable to human-designed networks for object recognition. Similarly, **LEAF (Learning Evolutionary AI Framework)** by Liang et al. (2019) [993] applied evolutionary algorithms to optimize both the architecture and hyperparameters of deep neural networks, demonstrating state-of-the-art performance in real-world tasks such as medical imaging and natural language processing. Meanwhile, Vargas and Murata (2016) [994] proposed **Spectrum-Diverse Neuroevolution**, which introduces a diversity-preserving mechanism to enhance the robustness of
evolved networks by maintaining variation at the behavioral level. These approaches highlight how evolutionary methods can lead to automated and efficient deep learning model discovery, reducing human intervention while improving performance.

Beyond architecture search, researchers have explored evolutionary algorithms as an alternative to gradient-based training. Such et al. (2017) [995] demonstrated that simple genetic algorithms could rival traditional methods like Q-learning and policy gradients in reinforcement learning tasks, particularly in environments with sparse rewards. This work challenged the conventional reliance on backpropagation, showing that evolutionary strategies could efficiently explore high-dimensional parameter spaces. More recently, **Fast-DENSER** by Assunção et al. (2021) [996] leveraged grammatical evolution to represent and evolve deep neural network structures, allowing for an efficient search over complex architectures with reduced computational overhead. Additionally, **Interactively Constrained Neuro-Evolution (ICONE)** by Rempis (2012) [997] introduced constraint masks to restrict the search space, enabling the evolution of highly specialized neural controllers while incorporating domain knowledge. These contributions collectively demonstrate the power of neuro-genetic evolution in optimizing deep learning architectures and training strategies, providing robust alternatives to traditional approaches.

<table><caption>Table 16.51: Summary of Contributions in Neuro-Genetic Evolution Research</caption><thead><tr><td>Authors (Year)</td><td>Contribution</td></tr></thead><tbody><tr><td>Stanley & Miikkulainen (2002)</td><td>Introduced NeuroEvolution of Augmenting Topologies (NEAT), a genetic algorithm that evolves both the structure and weights of neural networks. NEAT starts with minimal architectures and gradually complexifies them while maintaining diversity through speciation.</td></tr><tr><td>Stanley et al. (2005)</td><td>Extended the NEAT framework, demonstrating its ability to evolve increasingly complex and high-performing neural architectures through systematic augmentation.</td></tr><tr><td>Gauci & Stanley (2007)</td><td>Developed HyperNEAT, which leverages compositional pattern-producing networks (CPPNs) to encode connectivity patterns, enabling the evolution of large-scale networks with inherent symmetries and regularities.</td></tr><tr><td>Kassahun & Sommer (2005)</td><td>Proposed Evolutionary Acquisition of Neural Topologies (EANT), where neural networks start simple and gain complexity over time, reducing computational costs while improving network efficiency.</td></tr><tr><td>Miikkulainen et al. (2024)</td><td>Introduced CoDeepNEAT, an extension of NEAT that optimizes deep learning architectures by evolving topology, hyperparameters, and components, achieving performance comparable to human-designed networks.</td></tr><tr><td>Liang et al. (2019)</td><td>Developed the Learning Evolutionary AI Framework (LEAF), which applies evolutionary algorithms to optimize both neural architectures and hyperparameters, achieving state-of-the-art results in medical imaging and natural language processing.</td></tr><tr><td>Vargas & Murata (2016)</td><td>Proposed Spectrum-Diverse Neuroevolution, a method that preserves diversity at the behavioral level to enhance the robustness of evolved networks.</td></tr></tbody></table>

*Continued on next page*
(Continued from previous page)

<table><thead><tr><td>Authors (Year)</td><td>Contribution</td></tr></thead><tbody><tr><td>Such et al. (2017)</td><td>Demonstrated that simple genetic algorithms can rival traditional gradient-based approaches such as Q-learning and policy gradients in reinforcement learning tasks, particularly in sparse reward environments.</td></tr><tr><td>Assunção et al. (2021)</td><td>Developed Fast-DENSER, which utilizes grammatical evolution to efficiently search for deep neural network architectures while reducing computational overhead.</td></tr><tr><td>Rempis (2012)</td><td>Introduced Interactively Constrained Neuro-Evolution (ICONE), incorporating constraint masks to restrict the search space and evolve specialized neural controllers with domain-specific knowledge.</td></tr></tbody></table>

### 16.5.2.2 Recent Literature Review of Neuro-Genetic Evolution

Deep Learning Neuro-Genetic Evolution represents a significant fusion of evolutionary computation with deep learning, enhancing neural network optimization through genetic algorithms. Several key contributions have advanced this field by integrating adaptive genetic selection mechanisms into deep learning frameworks, allowing models to evolve dynamically. For instance, Stanley et. al. (2019) [998] pioneered NeuroEvolution of Augmenting Topologies (NEAT), which enables artificial neural networks to evolve both in structure and parameters, outperforming traditional gradient-based methods in complex problem-solving scenarios. Similarly, Bertens and Lee (2019) [999] examined the synergies between neural networks and evolutionary algorithms, emphasizing their application in biological neural modeling and robotics, showcasing how nature-inspired selection can enhance artificial intelligence adaptability. Expanding upon these foundations, Wang et al. (2023) [1000] demonstrated the efficacy of genetic evolution in feature selection for bioinformatics, where evolutionary deep learning significantly improved genomic data classification and disease prediction accuracy.

Another groundbreaking application of neuro-genetic learning is in autonomous systems and robotics, where deep reinforcement learning benefits from genetic optimization strategies. Pagliuca et. al. (2020) [1001] proposed a neuro-evolutionary approach that refines robotic movement through self-adaptive evolutionary deep learning, demonstrating superior efficiency in industrial automation and autonomous decision-making. Behjat et. al. (2019) [1002] introduced AGENT, a neuro-evolution framework that evolves both the topology and weights of neural networks. This adaptive approach addresses issues like premature convergence and stagnation, leading to improved performance in autonomous systems. The framework's efficacy is demonstrated through applications in unmanned aerial vehicle (UAV) collision avoidance, showcasing its relevance to industrial automation and autonomous decision-making. Similarly, Ahmed et. al. (2023) [1003] introduced a genetic reinforcement learning adaptation mechanism, optimizing Deep Q-networks for real-time applications such as self-driving cars and robotic control. The study illustrated how genetic selection enables reinforcement learning models to adapt faster to highly dynamic environments. Further, Miikkulainen et al. (2023) [1004] tackled the challenge of hyperparameter tuning in deep learning, leveraging evolutionary strategies to automate model selection and improve computational efficiency. Their work indicated that genetic-based hyperparameter optimization can outperform traditional grid and random search methods by at least 40 percentage in training speed and accuracy.

The impact of Neuro-Genetic Evolution extends into fields like cybersecurity, finance, and healthcare, where evolving models can optimize performance in unpredictable environments. Kannan et al. (2024) [1005] presented a neuro-genetic deep learning framework for IoT security, effectively detecting RPL attacks through an adaptive, self-learning anomaly detection system. This under-
scores the real-time adaptability of evolutionary deep learning in critical security applications. In financial forecasting, Zeng et. al. (2022) [1006] developed a hybrid genetic-deep learning model for predicting market fluctuations, demonstrating higher accuracy and robustness in volatile financial trends. Lastly, the use of Neuro-Genetic Evolution in software engineering and machine learning automation is gaining traction. S KV and Swamy (2024) [1007] explored ensemble-based neuro-genetic models to improve software quality by refining feature selection and defect prediction. By integrating genetic evolution strategies with deep learning ensembles, they significantly enhanced the reliability and performance of software engineering models. As deep learning systems become increasingly complex, leveraging evolutionary genetic strategies ensures their continual adaptation and optimization, pushing the boundaries of intelligent automation, real-time decision-making, and computational efficiency. These contributions collectively indicate that the fusion of deep learning with evolutionary computation is not only revolutionizing neural network architectures but also paving the way for next-generation AI models that can autonomously evolve and self-optimize across various domains.

<table><caption>Table 16.52: Summary of Recent Contributions in Neuro-Genetic Evolution Research</caption><thead><tr><td>Authors (Year)</td><td>Contribution</td></tr></thead><tbody><tr><td>Stanley et al. (2019)</td><td>Pioneered NeuroEvolution of Augmenting Topologies (NEAT), enabling artificial neural networks to evolve in structure and parameters, outperforming traditional gradient-based methods.</td></tr><tr><td>Bertens and Lee (2019)</td><td>Explored synergies between neural networks and evolutionary algorithms, particularly in biological neural modeling and robotics, showcasing adaptability improvements in AI through nature-inspired selection.</td></tr><tr><td>Wang et al. (2023)</td><td>Demonstrated genetic evolution in feature selection for bioinformatics, improving genomic data classification and disease prediction accuracy through evolutionary deep learning.</td></tr><tr><td>Pagliuca et al. (2020)</td><td>Proposed a neuro-evolutionary approach for robotic movement refinement via self-adaptive evolutionary deep learning, enhancing efficiency in industrial automation and autonomous decision-making.</td></tr><tr><td>Behjat et al. (2019)</td><td>Introduced AGENT, a neuro-evolution framework evolving both topology and weights of neural networks to prevent premature convergence and stagnation, applied to UAV collision avoidance.</td></tr><tr><td>Ahmed et al. (2023)</td><td>Developed a genetic reinforcement learning mechanism to optimize Deep Q-networks for real-time applications such as self-driving cars and robotic control.</td></tr><tr><td>Miikkulainen et al. (2023)</td><td>Applied evolutionary strategies for hyperparameter tuning in deep learning, achieving a 40% improvement in training speed and accuracy compared to traditional methods.</td></tr><tr><td>Kannan et al. (2024)</td><td>Designed a neuro-genetic deep learning framework for IoT security, enabling adaptive, self-learning anomaly detection systems for RPL attack detection.</td></tr><tr><td>Zeng et al. (2022)</td><td>Developed a hybrid genetic-deep learning model for financial forecasting, demonstrating increased accuracy and robustness in predicting volatile market trends.</td></tr></tbody></table>
<table><tr><td>S KV and Swamy (2024)</td><td>Explored ensemble-based neuro-genetic models for software quality improvement, refining feature selection and defect prediction using genetic evolution strategies.</td></tr></table>

### 16.5.2.3 Mathematical Analysis of Neuro-Genetic Evolution

The Neuro-Genetic Evolution Method in Artificial Intelligence (AI) integrates artificial neural networks (ANNs) and genetic algorithms (GAs) to optimize network structures, parameters, and learning strategies. This approach mimics biological evolution by iteratively refining a population of neural networks through genetic operations such as selection, crossover, and mutation. The

<figure><img src="image_212.png" /><figcaption>Figure 16.44: Visualization of Neuro-Genetic Evolution (NGE). Each panel shows a neural network individual within the evolving population. Nodes are color-coded by type (blue = input, green = hidden, red = output). Edges are weighted connections, with blue arrows representing positive weights and red arrows negative weights. Thickness corresponds to magnitude</figcaption></figure>

optimization process enhances both the architecture and the weights of the network, allowing it to learn complex patterns efficiently. Mathematically, let the neural network be represented as a function $f: \mathbb{R}^n \to \mathbb{R}^m$, where $n$ is the number of inputs, and $m$ is the number of outputs. The function is defined as:

$$ \mathbf{y} = f(\mathbf{x}, \mathbf{W}) = \sigma(\mathbf{W}\mathbf{x} + \mathbf{b}) \quad (16.360) $$

where $\mathbf{x} \in \mathbb{R}^n$ is the input vector, $\mathbf{W}$ is the weight matrix, $\mathbf{b}$ is the bias vector, and $\sigma(\cdot)$ is an activation function. The objective is to optimize $\mathbf{W}$ and $\mathbf{b}$ using genetic evolution. A population of neural networks is represented as:

$$ \mathcal{P} = \{(\mathbf{W}_1, \mathbf{b}_1), (\mathbf{W}_2, \mathbf{b}_2), \dots, (\mathbf{W}_N, \mathbf{b}_N)\} \quad (16.361) $$

where $N$ is the population size. The fitness function $F: \mathcal{P} \to \mathbb{R}$ evaluates the performance of each network on a given task:

$$ F(\mathbf{W}_i, \mathbf{b}_i) = - \sum_{j=1}^{M} \| f(\mathbf{x}_j, \mathbf{W}_i) - \mathbf{y}_j \|^2 \quad (16.362) $$

where $M$ is the number of training samples, $(\mathbf{x}_j, \mathbf{y}_j)$ are training pairs, and the objective is to maximize $F$. Selection is performed using a probabilistic function based on fitness:

$$ P_i = \frac{F(\mathbf{W}_i, \mathbf{b}_i)}{\sum_{k=1}^{N} F(\mathbf{W}_k, \mathbf{b}_k)} \quad (16.363) $$

where $P_i$ is the probability of selecting the $i$-th network for reproduction. Higher-fitness networks are more likely to be chosen. Crossover combines weights of two parent networks $(\mathbf{W}_A, \mathbf{b}_A)$ and $(\mathbf{W}_B, \mathbf{b}_B)$ to create an offspring $(\mathbf{W}_C, \mathbf{b}_C)$:

$$ \mathbf{W}_C = \alpha \mathbf{W}_A + (1 - \alpha) \mathbf{W}_B, \quad \mathbf{b}_C = \alpha \mathbf{b}_A + (1 - \alpha) \mathbf{b}_B \quad (16.364) $$
where $\alpha \sim U(0, 1)$ is a random crossover parameter. Mutation perturbs the weights and biases with a small random noise:

$$ \mathbf{W}' = \mathbf{W} + \eta \mathcal{N}(0, \sigma^2), \quad \mathbf{b}' = \mathbf{b} + \eta \mathcal{N}(0, \sigma^2) \qquad (16.365) $$

where $\eta$ is the mutation rate and $\mathcal{N}(0, \sigma^2)$ is Gaussian noise. The evolution process iterates over multiple generations $G$ to refine the population:

$$ \mathcal{P}^{(t+1)} = \text{Mutation}(\text{Crossover}(\text{Selection}(\mathcal{P}^{(t)}))) \qquad (16.366) $$

where $t$ denotes the generation index. Through this iterative process, neural networks evolve toward an optimal configuration, achieving better generalization and adaptation to their learning tasks. The convergence of the method is influenced by hyperparameters such as population size $N$, mutation rate $\eta$, and crossover probability $p_c$, which are optimized to balance **exploration** and **exploitation** in the evolutionary process.

### 16.5.2.4 Python Code to generate Figure 16.44

Neuro-Genetic Evolution (NGE) is usually represented by showing how a population of neural networks (genomes) evolves under genetic algorithms:

1. Each individual is a neural network (nodes + weighted connections).

2. A population of these networks exists at each generation.

3. Crossover and mutation create new networks.

4. Selection propagates the fittest.

To visualize NGE, we can plot:

1. Population at one generation → networks side by side.

2. Evolutionary flow → how parents produce offspring.

The Python code below produces the Figure 16.44 illustrating the Neuro-Genetic Evolution (NGE).

```python
import matplotlib.pyplot as plt
import networkx as nx
import random

# ----------------------------
# Create a random neural network
# ----------------------------
def build_random_network(nInputs=2, nHidden=2, nOutputs=1):
    G = nx.DiGraph()

    # Assign positions for visualization
    positions = {}
    idx = 0

    # Input layer
    for i in range(nInputs):
        node = f"I{i+1}``
        G.add_node(node, layer="input", pos=(i, 0))
        idx += 1

    # Hidden layer
```
```python
for i in range(n_hidden):
    node = f"H{i+1}"
    G.add_node(node, layer="hidden", pos=(i, 1))
    idx += 1

# Output layer
for i in range(n outputs):
    node = f"O{i+1}"
    G.add_node(node, layer="output", pos=(i, 2))
    idx += 1

# Random connections
nodes = list(G.nodes())
for _ in range(randomrandint(3, 6)):
    u, v = random.sample(nodes, 2)
    if u != v:
        G.add_edge(u, v, weight =random.random() uniform(-1, 1))

return G

# ------------------------
# Draw one network
# ------------------------
def draw_network(G, ax, title=""):
    pos = {n: d["pos"] for n, d in G.nodes(data=True)}

    # Node colors by layer
    colors = []
    for n, d in G.nodes(data=True):
        if d["layer"] == "input":
            colors.append("lightblue")
        elif d["layer"] == "hidden":
            colors.append("lightgreen")
        else:
            colors.append("salmon")

    nx.draw_networkx_nodes(G, pos, node_size=800, node_color=colors,
                           edgecolors="k", ax=ax)
    nx.draw_networkx labels(G, pos, ax=ax, font_size=8)

    # Draw weighted edges
    for u, v, d in GEdges(data=True):
        w = d["weight"]
        color = "blue" if w > 0 else "red"
        width = 1 + 2 * abs(w)
        ax annotate("",
                   xy =pos[v], xycoordinates="data",
                   xytext =pos[u], textcoordinates="data",
                   arrowprops =dict(arrowstyle="->",
                                     color =color,
                                     linewidth =width,
                                     alpha =0.6))
    ax.set_title(title, fontsize=9)
    ax .axis("off")
```
```python
# Plot NGE population
def plot_ngePopulation(pop_size=4, save_as=None):
    fig, axes = plt.subplots(1, pop_size, figsize=(4*pop_size, 3))
    for i in range(pop_size):
        G = build_random_network()
        draw_network(G, axes[i], title=f"Individual {i+1}")
    plt.suptitle("Neuro-Genetic Evolution: Population of Neural Networks", fontsize=14)
    if save_as:
        plt.savefig save_as, dπ=300, bbox_inches="tight")
    plt.show()

# Example run
if __name__ == "__main__":
    plot_ngePopulation(pop_size=4, save_as="nge population.png")
```

## 16.5.3 Cellular Encoding (CE)

### 16.5.3.1 Literature Review of Cellular Encoding (CE)

Cellular Encoding (CE) has played a significant role in the evolution of artificial neural networks by providing a structured and biologically inspired approach to encoding network architectures as labeled trees. The foundational work by Gruau (1993) [1008] introduced Cellular Encoding as a graph grammar-based method, which facilitates the application of genetic algorithms to evolve complex neural networks. This work was instrumental in demonstrating how CE can represent neural structures efficiently, enabling evolutionary optimization to generate compact and generalizable architectures. Following this, Gruau et. al. (1996) [1009] compared CE with direct encoding methods and provided compelling evidence that CE-based representations lead to superior neural architectures, particularly in terms of scalability and structural efficiency. This comparison underscored the advantage of hierarchical, modular encoding schemes over traditional direct encoding approaches, which often suffered from excessive parameterization and computational overhead. Furthermore, the application of CE to neurocontrol, as explored in Gruau and Whitley (1993) [1010], showcased its practical utility in dynamic system control, emphasizing the potential of evolved neural controllers in real-world tasks requiring adaptive and self-organizing behavior.

A critical examination of the structural capacity of CE was conducted by Gutierrez et. al. (2004) [1011], who investigated its ability to generate a diverse range of feedforward neural network topologies. Their study provided empirical validation of CE's robustness in evolving networks with varying complexities, highlighting its versatility in neural architecture search. Meanwhile, Zhang and Muhlenbein (1993) [1012] examined the role of Occam's Razor in CE-based neural evolution, revealing how CE can balance model complexity and performance, leading to optimized network topologies that avoid unnecessary overfitting. This work reinforced the idea that evolutionary principles, when combined with structured encoding mechanisms, can yield networks that are both minimalistic and effective. On a related note, Kitano (1990) [1013] introduced a graph generation system for designing neural networks using genetic algorithms, which, while predating formal CE methodologies, provided foundational insights into evolutionary network design principles. Kitano's work laid the groundwork for later advancements in CE by demonstrating the feasibility of using genetic repre-
sentations to encode and evolve neural architectures efficiently.

A significant parallel development came with Miller and Turner (2015) [1014] and Miller (2020) [1015]. Their work on Cartesian Genetic Programming (CGP), which, although not directly related to CE, introduced an alternative method of encoding neural architectures using graph-based representations. CGP's emphasis on modular and reusable computational structures provided complementary insights into encoding strategies for neural networks, influencing later refinements in CE-based models. Similarly, Stanley and Miikkulainen (2002) [983] introduced the NEAT (NeuroEvolution of Augmenting Topologies) algorithm, which evolved neural networks by progressively augmenting topologies. While NEAT follows a different evolutionary strategy than CE, it shares the core philosophy of dynamically evolving neural structures rather than relying on fixed architectures. The interplay between NEAT and CE highlights the broader landscape of evolutionary neural architecture search, where different encoding schemes lead to varying trade-offs in efficiency, scalability, and adaptability.

In more recent advancements, Hernandez Ruiz et al. (2021) [1016] extended the CE paradigm to neural cellular automata (NCA), proposing a manifold representation capable of generating diverse images through dynamic convolutional mechanisms. This work demonstrated how CE-inspired principles could be applied beyond classical neural network design, expanding its applicability to generative models and self-organizing systems. Furthermore, Hajij, Istvan, and Zamzmi (2020) [1017] introduced Cell Complex Neural Networks (CXNs), which generalized message-passing schemes to higher-dimensional structures, offering a novel perspective on encoding neural computations. Their approach emphasized the mathematical rigor behind encoding mechanisms and provided new avenues for incorporating CE-like principles into modern deep learning architectures.

Overall, these works collectively underscore the versatility and depth of Cellular Encoding in neural network evolution. From its inception as a biologically motivated encoding scheme to its applications in neurocontrol, architecture search, and complex generative models, CE has demonstrated its capability to produce efficient, scalable, and structurally diverse networks. The cross-pollination of ideas between CE, CGP, NEAT, and other evolutionary strategies further enriches the field, suggesting that hybrid approaches leveraging the strengths of multiple encoding mechanisms may pave the way for the next generation of evolutionary neural networks. Through rigorous analysis and empirical validation, these studies illustrate how structured encoding methods can significantly enhance neural network optimization, making CE a foundational pillar in the ongoing advancement of neuroevolution.

<table><caption>Table 16.53: Summary of Contributions in Cellular Encoding Research</caption><thead><tr><td>Authors (Year)</td><td>Contribution</td></tr></thead><tbody><tr><td>Gruau (1993)</td><td>Introduced Cellular Encoding (CE) as a graph grammar-based approach to representing neural networks. Demonstrated how CE can facilitate the evolution of complex neural structures through genetic algorithms.</td></tr><tr><td>Gruau (1996)</td><td>Compared Cellular Encoding with direct encoding methods, proving CE's superiority in producing compact, efficient, and generalizable neural networks. Showed how hierarchical representations improve evolutionary search.</td></tr><tr><td>Gruau & Whitley (1993)</td><td>Applied Cellular Encoding to neurocontrol problems, illustrating its ability to evolve adaptive and robust neural controllers for dynamic systems.</td></tr></tbody><tfoot><tr><td colspan="2">Continued on next page</td></tr></tfoot></table>
(Continued from previous page)

<table><thead><tr><td>Authors (Year)</td><td>Contribution</td></tr></thead><tbody><tr><td>Gutierrez et. al. (2004)</td><td>Investigated the capacity of CE to generate diverse feedforward neural network topologies, providing empirical evidence of its flexibility and scalability.</td></tr><tr><td>Zhang & Muhlenbein (1993)</td><td>Explored the application of CE with genetic algorithms under Occam's Razor, demonstrating how CE can evolve minimal yet high-performing neural architectures.</td></tr><tr><td>Kitano (1990)</td><td>Introduced a graph generation system for evolving neural networks using genetic algorithms, laying a theoretical foundation for Cellular Encoding approaches.</td></tr><tr><td>Miller & Turner (2015) and Miller (2020)</td><td>Developed Cartesian Genetic Programming (CGP), a related graph-based encoding approach that influenced CE by reinforcing modular and reusable neural components.</td></tr><tr><td>Stanley & Miikku-lainen (2002)</td><td>Proposed the NEAT algorithm, which evolves neural networks through progressive augmentation of topologies. Though distinct from CE, NEAT shares principles of evolving efficient and adaptable architectures.</td></tr><tr><td>Hernandez Ruiz et al. (2021)</td><td>Extended Cellular Encoding to Neural Cellular Automata (NCA), demonstrating its application in generative models and image synthesis using convolutional mechanisms.</td></tr><tr><td>Hajij, Istvan, & Zamzmi (2020)</td><td>Introduced Cell Complex Neural Networks (CXNs), generalizing message-passing schemes to higher-dimensional structures, providing a new theoretical framework for encoding neural computations.</td></tr></tbody></table>

### 16.5.3.2 Recent Literature Review of Cellular Encoding (CE)

Cellular Encoding (CE) in neural networks has emerged as a pivotal concept bridging computational neuroscience, bioinformatics, and artificial intelligence. At its core, CE leverages biologically inspired encoding strategies to represent, manipulate, and interpret cellular processes and neural dynamics. Sun et al. (2025) [1018] explored this phenomenon by investigating how learning transforms hippocampal neural representations into an orthogonalized state machine. Their study found that as learning progresses, cellular and population-level neural activity becomes increasingly structured, forming distinct encoded representations that optimize memory retrieval and spatial navigation. By linking artificial neural networks (ANNs) with hippocampal activity, their work provides a foundational framework for biologically plausible machine learning models. Similarly, Hu et al. (2025) [878] extended this idea to genomics by designing an ensemble deep learning framework for long non-coding RNA (lncRNA) subcellular localization, demonstrating how CE can aid in deciphering complex gene regulation patterns.

Advancing the scope of CE, Guan et al. (2025) [1019] developed a graph neural structure encoding approach for semantic segmentation of nuclei in pathological tissues. Their work leverages graph neural networks (GNNs) to model spatial relationships between cellular structures, allowing precise delineation of subcellular components. This marks a critical step in biomedical imaging, where accurate segmentation is essential for disease diagnosis and histopathological analysis. On the computational front, Ghosh et al. (2025) [1020] tackled regulatory network encoding by designing a deep learning framework for transcription factor binding site prediction. Their work integrated DNABERT and convolutional neural networks (CNNs) to extract and encode DNA sequence motifs, significantly enhancing the accuracy of binding site identification. This study exemplifies how CE techniques can revolutionize functional genomics, enabling more efficient identification of ge-
netic regulatory elements.

Beyond genomics, CE also finds applications in cellular perturbation modeling and neuroinformatics. Sun et al. (2025) [1021] introduced a perturbation proteomics-based virtual cell model that encodes dynamic protein interaction networks. By integrating large-scale perturbation data with deep learning architectures, their approach provides a novel way to simulate cellular responses to environmental and pharmacological interventions. Grosjean et al. (2025) [1022] contributed to self-supervised learning in neuroscience by developing a network-aware encoding strategy for detecting genetic modifiers of neuronal activity. Their work highlights the importance of high-content phenotypic screening, where CE can be used to identify how genetic variations influence neural dynamics. These studies collectively illustrate CE's potential in predictive modeling of cellular behaviors, from molecular interactions to large-scale neuronal networks.

Expanding CE applications to neurodevelopmental disorders, de Carvalho et al. (2025) [1024] conducted a gene network analysis of autism spectrum disorder (ASD) by encoding synaptic and cellular alterations. Their work linked transcription factor mutations to widespread disruptions in neuronal communication, paving the way for targeted therapeutic strategies. Similarly, Gonzalez et al. (2025) [1023] introduced an in vivo single-cell electroporation tool that enables real-time encoding of hippocampal neurons. By integrating genetically encoded calcium indicators and voltage sensors, their method provides unprecedented control over synaptic plasticity and neural circuit modulation. These innovations reflect the expanding role of CE in understanding and manipulating brain function at a cellular level.

Lastly, broader implications of CE extend to brain connectivity and non-neuronal function. Sprecher (2025) [1025] investigated how neural networks encode and coordinate brain-wide activity, providing a computational model of synaptic connectivity. Their findings emphasize the role of CE in shaping neural excitability and dynamic regulation. Li et al. (2025) [1026] expanded this perspective by examining non-neuronal contributions to neural encoding. Their work revealed that glial cells and extracellular matrix components actively participate in shaping encoded neural signals, challenging the long-held neuron-centric view of the nervous system. These studies reinforce the fundamental principle that CE is not confined to neurons alone but encompasses an intricate interplay of cellular components, making it a cornerstone of modern neuroscience and bioinformatics.

Together, these studies form a comprehensive landscape of Cellular Encoding, demonstrating its versatility in biological data representation, computational modeling, and disease research. Whether through neural network-based segmentation, regulatory sequence analysis, or predictive perturbation modeling, CE is driving innovation across multiple disciplines. The convergence of deep learning, systems biology, and neuroscience highlights CE's ability to bridge biological complexity with computational efficiency, ultimately leading to smarter AI models, better disease diagnostics, and deeper insights into cellular intelligence.

<table><caption>Table 16.54: Summary of Recent Contributions in Cellular Encoding Research</caption><thead><tr><td>Authors (Year)</td><td>Contribution</td></tr></thead><tbody><tr><td>Sun et al. (2025)</td><td>Investigated how learning transforms hippocampal neural activity into an orthogonalized state machine, demonstrating how structured encoding optimizes memory retrieval and spatial navigation.</td></tr><tr><td>Hu et al. (2025)</td><td>Developed an ensemble deep learning framework for long non-coding RNA (lncRNA) subcellular localization, showcasing how cellular encoding enhances gene regulation analysis.</td></tr></tbody></table>

Continued on next page
Continued from previous page

<table><thead><tr><td>Authors (Year)</td><td>Contribution</td></tr></thead><tbody><tr><td>Guan et al. (2025)</td><td>Proposed a graph neural structure encoding method for semantic segmentation of nuclei in pathological tissues, enabling improved cellular structure identification in biomedical imaging.</td></tr><tr><td>Ghosh et al. (2025)</td><td>Designed a deep learning-based transcription factor binding site predictor using DNABERT and convolutional neural networks (CNNs) to extract and encode DNA sequence motifs.</td></tr><tr><td>Sun et al. (2025)</td><td>Introduced a perturbation proteomics-based virtual cell model, integrating protein interaction networks with deep learning to simulate cellular responses under environmental and pharmacological perturbations.</td></tr><tr><td>Grosjean et al. (2025)</td><td>Developed a self-supervised learning approach for detecting genetic modifiers of neuronal activity, using a network-aware encoding strategy to enhance high-content phenotypic screening.</td></tr><tr><td>de Carvalho et al. (2025)</td><td>Conducted a gene network analysis on autism spectrum disorder (ASD), encoding synaptic and cellular alterations linked to transcription factor mutations affecting neuronal communication.</td></tr><tr><td>Gonzalez et al. (2025)</td><td>Created an in vivo single-cell electroporation tool, enabling real-time encoding of hippocampal neurons through genetically encoded calcium indicators and voltage sensors.</td></tr><tr><td>Sprecher (2025)</td><td>Investigated how neural networks encode and regulate brain-wide connectivity, providing a computational model for synaptic excitability and neural dynamics.</td></tr><tr><td>Li et al. (2025)</td><td>Examined non-neuronal contributions to neural encoding, revealing that glial cells and extracellular matrix components play an active role in shaping encoded neural signals.</td></tr></tbody></table>

### 16.5.3.3 Mathematical Analysis of Cellular Encoding (CE)

Cellular Encoding (CE) is a powerful and biologically inspired method for encoding artificial neural networks (ANNs) in evolutionary computation. It is based on the principle that a genotype (a compact representation) can be mapped to a phenotype (a fully structured ANN) through a set of developmental rules. These rules mimic the way biological organisms develop from a single cell through cellular divisions and differentiations governed by genetic instructions. The CE method introduces a recursive process, where the construction of an ANN follows a sequence of genetic instructions that define cell divisions, differentiations, and synaptic connections.

Mathematically, let $G$ be the genotype, which consists of a set of developmental instructions $\{g_i\}$, where each $g_i$ is a rule that governs a cellular transformation. The phenotype $P$, which is the neural network, is generated through a function $\mathcal{D}$ that applies these instructions iteratively:

$$P = \mathcal{D}(G) \qquad (16.367)$$

where $\mathcal{D}$ is a mapping function that executes a sequence of operations to form the final neural network. Each cell in this developmental process can be represented by a state $S$, which consists of attributes such as its position $x$, its type $t$, and its connectivity $C$:

$$S_i = (x_i, t_i, C_i) \qquad (16.368)$$

where $x_i \in \mathbb{R}^n$ represents the spatial coordinates of the cell, $t_i$ denotes the type of the neuron (such as input, hidden, or output), and $C_i$ represents the set of synaptic connections formed during
development. The developmental process begins with a single root cell at an initial state $S_0$. The recursive application of genetic instructions modifies this state according to transformation rules. Each transformation is mathematically represented by a function $T$, such that at step $k$:

$$S_i^{(k+1)} = T(S_i^{(k)}, g_k) \qquad (16.369)$$

where $g_k$ specifies an operation such as cell division, differentiation, or connection formation. Cell division is a key operation in CE and can be modeled as a function that creates two daughter cells $S_A$ and $S_B$ from a parent cell $S_P$:

$$S_A, S_B = \Phi(S_P, g_d) \qquad (16.370)$$

where $g_d$ is a division instruction. The transformation function $\Phi$ ensures that spatial attributes

Cellular Encoding (CE) Derivation Tree

<figure><img src="image_213.png" /><figcaption>Figure 16.45: Cellular Encoding (CE) derivation tree. A developmental grammar tree specifies transformations</figcaption></figure>

$x$ and connectivity $C$ are updated to reflect the new structure:

$$x_A = x_P + \delta x, \quad x_B = x_P - \delta x \qquad (16.371)$$

$$C_A = f(C_P), \quad C_B = g(C_P) \qquad (16.372)$$

where $\delta x$ is a spatial displacement vector, and $f, g$ define how connectivity properties are inherited. Cell differentiation is governed by another transformation function $\Psi$, which modifies the type attribute:

$$t_A = \Psi(t_P, g_t) \qquad (16.373)$$

where $g_t$ is a differentiation rule that determines whether a neuron becomes an excitatory, inhibitory, or output neuron. Synaptic connectivity is established through a connection function $\Gamma$, which creates weighted edges between neurons:

$$C_{ij} = \Gamma(S_i, S_j, g_c) \qquad (16.374)$$
<figure><img src="image_214.png" /><figcaption>Figure 16.46: Neural network produced by Cellular Encoding (CE) rules: The final evolved neural network with inputs (blue), hidden units (green), and outputs (red). Edge color encodes sign (blue = positive, red = negative), thickness encodes magnitude</figcaption></figure>

where $g_c$ is a connection rule that specifies how weights $w_{ij}$ are assigned based on distance and developmental constraints:

$$w_{ij} = \gamma \cdot e^{-\alpha \|x_i - x_j\|} \qquad (16.375)$$

where $\gamma$ is a scaling factor and $\alpha$ controls the decay of connection strength with distance. As the network grows, a hierarchical and structured topology emerges, leading to a functional ANN. The final output network is expressed as a weighted graph $\mathcal{N} = (V, E, W)$, where $V$ is the set of neurons, $E$ is the set of synaptic connections, and $W$ represents the weight matrix:

$$W = [w_{ij}]_{n \times n} \qquad (16.376)$$

The learning process in CE can also involve adaptive modifications, where the weights evolve according to a rule based on Hebbian learning:

$$\Delta w_{ij} = \eta \cdot x_i x_j \qquad (16.377)$$

where $\eta$ is the learning rate, ensuring that connectivity patterns are refined through evolutionary pressure. Additionally, mutation and crossover operators in genetic algorithms can modify $G$, leading to different developmental trajectories:

$$G' = \mathcal{M}(G), \quad G'' = \mathcal{C}(G_1, G_2) \qquad (16.378)$$

where $\mathcal{M}$ represents mutation and $\mathcal{C}$ represents crossover. Thus, the CE method provides a compact, recursive, and biologically inspired way to encode complex neural structures, ensuring efficient exploration of the search space for ANN architectures.
### 16.5.3.4 Python Code to generate Figures 16.45 and 16.46

Cellular Encoding (CE) is a fascinating neuroevolution approach. Instead of directly encoding weights like GNARL, CE uses a developmental grammar: a sequence of rules (a tree-like structure) that “grows” a neural network. To visualize it, we typically show:

1. A derivation tree (the CE program).

2. The neural network produced after applying the developmental rules.

The Python code below produces the Figures 16.45 and 16.46 illustrating the Cellular Encoding (CE).

```python
import matplotlib.pyplot as plt
import networkx as nx

# Example CE derivation tree

def build.ce_tree():
    """
    Example CE tree structure:
        - Start: Module
        - Split: produces hidden nodes
        - AddConn: adds connections
        - Recur: adds recurrent link
    """
    tree = nx.DiGraph()

    # Nodes in the CE program
    tree.add_node("Module", level=0)
    tree.add_node("Split", level=1)
    tree.add_node("AddConn", level=2)
    tree.add_node("Recur", level=2)

    # Edges of the derivation
    tree.add_edge("Module", "Split")
    tree.add_edge("Split", "AddConn")
    tree.add_edge("Split", "Recur")

    return tree

# Derived neural network

def build.ce_network():
    G = nx.DiGraph()

    # Input, hidden, output
    G.add_node("I1", layer="input", pos=(0,0))
    G.add_node("I2", layer="input", pos=(1,0))

    G.add_node("H1", layer="hidden", pos=(0.5,1))
    G.add_node("H2", layer="hidden", pos=(1.5,1))

    G.add_node("O1", layer="output", pos=(1,2))
```
```python
# Connections from "developmental rules"
G.add_edge("I1", "H1", weight=0.8)
G.add_edge("I2", "H2", weight=-0.6)
G.add_edge("H1", "O1", weight=0.9)
G.add_edge("H2", "O1", weight=0.5)
G.add_edge("H2", "H1", weight=0.4) # recurrent / lateral

return G
```
<figure><img src="image_215.png" /></figure>

16.5.4 GeNeralized Acquisition of Recurrent Links (GNARL)

16.5.4.1 Literature Review of GeNeralized Acquisition of Recurrent Links (GNARL)

The Generalized Acquisition of Recurrent Links (GNARL) algorithm, introduced by Saunders, Angeline, and Pollack (1993) [1027] represents a pioneering approach to the simultaneous evolution of both the topology and weights of recurrent neural networks (RNNs). Unlike traditional neural network training methods that rely on gradient-based optimization techniques such as backpropagation through time, GNARL employs an evolutionary strategy that allows for the dynamic adaptation of network structures. This early work established a robust foundation for the evolutionary acquisition of network architectures, demonstrating how an evolutionary algorithm could generate highly flexible, problem-specific RNN topologies. The authors showed that the algorithm could evolve networks with complex internal dynamics, crucial for tasks requiring memory and temporal dependencies. Expanding on these ideas, Angeline, Saunders, and Pollack (1994) [1028] further refined the approach by presenting an extensive comparison between GNARL and traditional methods, emphasizing the limitations of genetic algorithms in network evolution and advocating for the effectiveness of GNARL's approach in evolving RNNs more efficiently.

In subsequent years, research on neuroevolution continued to build on the principles introduced by GNARL, with Stanley and Miikkulainen (2002) [983] developing the NeuroEvolution of Augmenting Topologies (NEAT) framework, which introduced innovative mechanisms for evolving network topology while preserving functional structures. While NEAT diverged from GNARL in several ways—particularly in its use of speciation and historical markings—it acknowledged GNARL as one of the earliest attempts to evolve both the architecture and weights of RNNs. Similarly, Schmidhuber (1996) [1029] incorporated GNARL into his broader discussion on self-improving artificial intelligence systems, recognizing its capacity for autonomous structural adaptation. His work positioned GNARL within the broader context of incremental self-improvement methodologies, highlighting its potential for developing increasingly complex behaviors in multi-agent systems.

A critical review by Yao (1999) [1030] provided a comprehensive survey of artificial neural network evolution, situating GNARL within the broader landscape of evolutionary algorithms for neural networks. His work reinforced the significance of GNARL in pioneering structural evolution and weight adaptation, arguing that evolutionary approaches offered a promising alternative to traditional training methods, particularly for non-differentiable and highly non-linear optimization
problems. Further expanding on GNARL's contributions, Floreano, Dürr, and Mattiussi (2008) [1031] examined its role in the historical development of neuroevolution, emphasizing how its structural adaptability enabled more efficient learning processes compared to fixed-topology networks. Their work underlined GNARL's influence on later neuroevolutionary techniques that sought to balance exploration and exploitation in evolving network structures.

GNARL's applicability extended beyond theoretical discussions, influencing practical applications in reinforcement learning and control tasks. Gomez and Miikkulainen (1999) [1032] leveraged the principles of evolving recurrent neural networks, as demonstrated by GNARL, to solve non-Markovian control tasks using neuroevolution. Their work underscored the importance of evolving memory-capable networks, an area where GNARL had previously demonstrated efficacy. Similarly, Kassahun and Sommer (2005) [987] acknowledged GNARL in their research on reinforcement learning, specifically in the context of evolutionary methods for optimizing neural network topologies. Their study demonstrated that evolutionary approaches, inspired by GNARL, could significantly improve the learning efficiency of reinforcement learning agents by allowing the network structure itself to adapt alongside the weight parameters.

The broader implications of GNARL's methodology can also be seen in Moriarty and Miikkulainen's (1996) [1033] exploration of reinforcement learning through symbiotic evolution. Their study built upon GNARL's principles by demonstrating how co-evolutionary strategies could lead to more efficient learning processes, particularly for tasks requiring coordination among multiple network components. Furthermore, Gomez and Miikkulainen (1997) [1034] extended these ideas by investigating the incremental evolution of complex behaviors, a concept deeply rooted in GNARL's evolutionary framework. Their research emphasized the importance of evolving modular and hierarchical structures, recognizing GNARL's role in shaping later work on adaptive network evolution. Collectively, these studies illustrate the enduring impact of GNARL on the field of neuroevolution, highlighting its foundational contributions to evolving RNNs for complex learning and control tasks.

<table><caption>Table 16.55: Summary of Contributions to GNARL</caption><thead><tr><td>Authors (Year)</td><td>Contribution</td></tr></thead><tbody><tr><td>Saunders, Angeline, and Pollack (1993)</td><td>Introduced GNARL as an evolutionary algorithm for evolving both the topology and weights of recurrent neural networks (RNNs). Demonstrated its ability to evolve unconstrained architectures with complex internal dynamics.</td></tr><tr><td>Angeline, Saunders, and Pollack (1994)</td><td>Further refined GNARL by comparing its efficiency with traditional methods such as genetic algorithms and gradient-based approaches. Highlighted its advantages in evolving structurally adaptive networks.</td></tr><tr><td>Stanley and Miikkulainen (2002)</td><td>Developed NEAT, an evolutionary algorithm for augmenting neural network topologies. Recognized GNARL as an early attempt at evolving both structure and weights, influencing later neuroevolution research.</td></tr><tr><td>Schmidhuber (1996)</td><td>Discussed GNARL in the context of self-improving AI systems. Positioned GNARL within the framework of incremental self-improvement and multi-agent learning.</td></tr></tbody></table>

*Continued on next page*
*Continued from previous page*

<table><thead><tr><td>Authors (Year)</td><td>Contribution</td></tr></thead><tbody><tr><td>Yao (1999)</td><td>Provided a comprehensive review of evolutionary artificial neural networks, citing GNARL as a foundational approach for evolving both network structure and parameters. Highlighted its applicability to non-differentiable optimization problems.</td></tr><tr><td>Floreano, Dürr, and Matiussi (2008)</td><td>Examined GNARL's role in the historical development of neuroevolutionary methods. Emphasized its ability to balance exploration and exploitation in evolving network structures.</td></tr><tr><td>Gomez and Miikkulainen (1999)</td><td>Applied GNARL-inspired neuroevolution techniques to solve non-Markovian control tasks. Showed the importance of evolving memory-capable recurrent networks for reinforcement learning.</td></tr><tr><td>Kassahun and Sommer (2005)</td><td>Investigated reinforcement learning through evolutionary neural network optimization. Referenced GNARL as a key precursor to modern neuroevolutionary strategies.</td></tr><tr><td>Moriarty and Miikkulainen (1996)</td><td>Explored reinforcement learning using symbiotic evolution. Built on GNARL's principles to demonstrate how co-evolutionary strategies enhance learning in neural networks.</td></tr><tr><td>Gomez and Miikkulainen (1997)</td><td>Extended GNARL-based ideas to incremental evolution of complex behaviors. Emphasized the importance of evolving modular and hierarchical structures for adaptive learning.</td></tr></tbody></table>

### 16.5.4.2 Mathematical Analysis of Generalized Acquisition of Recurrent Links (GNARL)

The Generalized Acquisition of Recurrent Links (GNARL) method is a computationally advanced framework in artificial intelligence, specifically designed to evolve artificial neural networks (ANNs) dynamically. Unlike traditional gradient-based learning methods, GNARL employs evolutionary algorithms to optimize the network structure and weights, allowing for adaptive learning in recurrent neural networks (RNNs). The core principle behind GNARL is the use of genetic algorithms (GAs) to iteratively refine network topologies by introducing and modifying recurrent links based on a well-defined fitness function. Mathematically, this process can be understood through an iterative optimization problem where both the network connectivity *C* and weight matrices *W* evolve over time.

To represent an artificial neural network undergoing GNARL-based evolution, let *N* be the number of neurons in the network. The connectivity of the network at any given time step *t* is given by an adjacency matrix *C*(*t*), where each entry *C*<sub>ij</sub>(*t*) is defined as

$$C_{ij}(t) = \begin{cases} 1, & \text{if a connection exists from neuron } j \text{ to neuron } i \text{ at time } t, \\ 0, & \text{otherwise.} \end{cases} \qquad (16.379)$$

Each connection has an associated weight matrix *W*(*t*), whose elements *W*<sub>ij</sub>(*t*) evolve based on mutation and crossover mechanisms, following a genetic algorithm paradigm. The update rule for weights follows

$$W_{ij}(t+1) = W_{ij}(t) + \Delta W_{ij}(t), \qquad (16.380)$$
where the weight perturbation $\Delta W_{ij}(t)$ is given by

$$ \Delta W_{ij}(t) = \eta \frac{\partial F}{\partial W_{ij}} + \xi_{ij}, \qquad (16.381) $$

where $\eta$ is a learning rate, $F$ is a fitness function evaluating network performance, and $\xi_{ij}$ is a random mutation term typically drawn from a zero-mean Gaussian distribution $N(0, \sigma^2)$ to introduce stochasticity in the evolution process. The recurrent dynamics of the network are governed by the

<figure><img src="image_216.png" /><figcaption>Figure 16.47: Visualization of a GNARL (GeNeralized Acquisition of Recurrent Links) genome: Input nodes (blue), hidden nodes (green), and output nodes (red) are shown. Solid arrows represent standard feedforward connections, while dotted arrows indicate recurrent links. Edge color represents sign (blue: positive, red: negative), and thickness is proportional to weight magnitude</figcaption></figure>

activation function $\phi$, leading to neuron state updates given by

$$ x_i(t + 1) = \phi \left( \sum_{j=1}^{N} C_{ij}(t) W_{ij}(t) x_j(t) + b_i(t) \right), \qquad (16.382) $$

where $x_i(t)$ is the activation of neuron $i$ at time $t$, and $b_i(t)$ represents the neuron bias, which can also be evolved over time via

$$ b_i(t + 1) = b_i(t) + \Delta b_i(t), \qquad (16.383) $$

where

$$ \Delta b_i(t) = \eta \frac{\partial F}{\partial b_i} + \xi_i. \qquad (16.384) $$

The evolution of connectivity is achieved by a probabilistic mechanism where links are added or removed based on fitness evaluations. If $P_{add}$ and $P_{del}$ are the probabilities of adding and deleting
a link, respectively, then

$$P_{\text{add}} = \frac{1}{1 + e^{-\alpha(F-F_{\text{thresh}})}}, \quad P_{\text{del}} = 1 - P_{\text{add}}, \qquad (16.385)$$

where $\alpha$ is a sensitivity parameter and $F_{\text{thresh}}$ is a threshold fitness value that governs network complexity. The addition of a new connection follows

$$C_{ij}(t+1) = C_{ij}(t) + \mathbb{I}(U < P_{\text{add}}), \qquad (16.386)$$

where $U \sim \mathcal{U}(0, 1)$ is a uniformly distributed random variable and $\mathbb{I}$ is the indicator function. Similarly, link deletion follows

$$C_{ij}(t+1) = C_{ij}(t) - \mathbb{I}(U < P_{\text{del}}). \qquad (16.387)$$

The fitness function $F$ is domain-specific and often includes terms for error minimization, stability, and computational efficiency. A common formulation in supervised learning scenarios is

$$F = -\sum_{k=1}^{M} (y_k - \hat{y}_k)^2 + \lambda \sum_{i,j} C_{ij}(t), \qquad (16.388)$$

where $M$ is the number of training samples, $y_k$ and $\hat{y}_k$ are the target and predicted outputs, respectively, and $\lambda$ is a regularization parameter that penalizes excessive network complexity. The evolutionary process iterates through selection, mutation, and crossover operations. Selection follows a fitness-proportional rule given by

$$P_i = \frac{e^{\beta F_i}}{\sum_j e^{\beta F_j}}, \qquad (16.389)$$

where $\beta$ is a selection pressure parameter controlling the preference for higher-fitness individuals. Mutation is applied independently to weights and biases via

$$W_{ij}(t+1) = W_{ij}(t) + \mathcal{N}(0, \sigma_w^2), \quad b_i(t+1) = b_i(t) + \mathcal{N}(0, \sigma_b^2). \qquad (16.390)$$

Crossover occurs between two parent networks $A$ and $B$ to produce offspring $O$, where

$$W_{ij}^{O} = \begin{cases} W_{ij}^{A}, & \text{if } U < 0.5, \\ W_{ij}^{B}, & \text{otherwise} \end{cases} \qquad (16.391)$$

Recurrent links play a crucial role in GNARL's ability to discover temporal dependencies in sequential data. The recurrent connections enable information to persist across time steps, making the method particularly suited for problems such as time-series prediction and reinforcement learning. The recurrent update equations are formulated as

$$h_i(t+1) = \phi \left( \sum_j C_{ij}^{(R)}(t) W_{ij}^{(R)}(t) h_j(t) + \sum_j C_{ij}^{(I)}(t) W_{ij}^{(I)}(t) x_j(t) \right), \qquad (16.392)$$

where $C_{ij}^{(R)}$ and $W_{ij}^{(R)}$ denote recurrent connectivity and weights, while $C_{ij}^{(I)}$ and $W_{ij}^{(I)}$ denote input-to-hidden connections. The stability of the evolved networks is analyzed using Lyapunov exponents, where local stability requires

$$\max |\lambda_i(J(W, C))| < 1, \qquad (16.393)$$

where $J(W, C)$ is the Jacobian matrix of the system.
### 16.5.4.3 Python Code to generate Figure 16.47

GNARL is one of the earliest neuroevolution methods, focusing on evolving connection weights and recurrent links in neural networks. A visualization should emphasize:

1. Input, hidden, and output layers.

2. Recurrent links (self-loops or feedback connections).

3. Weighted edges (blue = positive, red = negative, width proportional to magnitude).

The Python code below produces the Figure 16.47 illustrating the GeNeralized Acquisition of Recurrent Links (GNARL).

```python
import matplotlib.pyplot as plt
import networkx as nx
import random

# Example GNARL genome

def generate_gnarl_genome(numInputs=3, numHidden=4, numOutputs=2,
                             seed=42):
    random.seed(seed)
    G = nx.DiGraph()

    # Input nodes
    inputs = [f"I{i}" for i in range(1, numInputs + 1)]
    for i, n in enumerate(Inputs):
        G.add_node(n, layer="input", pos=(i, 0))

    # Hidden nodes
    hidden = [f"H{i}" for i in range(1, numHidden + 1)]
    for i, n in enumerate(hidden):
        x, y = randomhea
        G.add_node(n, layer="hidden", pos=(x, y))

    # Output nodes
    outputs = [f"O{i}" for i in range(1, numOutputs + 1)]
    for i, n in enumerate(Outputs):
        G.add_node(n, layer="output", pos=(i, 2))

    # Input to hidden
    for inp in inputs:
        for hid in hidden:
            if random.random() < 0.7:
                G.add_edge inp, hid, weight=round(randomhea
2))

    # Hidden to output
    for hid in hidden:
        for out in outputs:
            if random.random() < 0.7:
                G.add_edge(hid, out, weight=round(randomhea
2))
```
```python
# Recurrent links (GNARL key feature)
for hid in hidden:
    if random.random() < 0.4:
        G.add_edge(hid, hid, weight=round(randomhea(-1, 1), 2))
        # self-loop
    for other in hidden:
        if hid != other and random.random() < 0.3:
            G.add_edge(hid, other, weight=round(randomhea(-1,
61
62
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

# ---------------------------
# Plotting function
# ---------------------------

def plot_gnarl(G, title="GNARL Genome", save_as=None):
    pos = {n: d["pos"] for n, d in G.nodes(data=True)}

    # Node colors
    color_map = []
    for n, d in G.nodes(data=True):
        if d["layer"] == "input":
            color_map.append("lightblue")
        elif d["layer"] == "hidden":
            color_map.append("lightgreen")
        else:
            color_map.append("salmon")

    # Draw nodes
    nx.draw_networkx_nodes(G, pos, node_color=color_map, node_size=800,
                           edgecolors="k")
    nx.draw_networkx labels(G, pos, font_size=8)

    # Draw edges (with weight-based styling)
    ax = plt.gca()
    for u, v, d in GEdges(data=True):
        w = d["weight"]
        color = "blue" if w > 0 else "red"
        width = 1 + 2 * abs(w)
        alpha = 0.7
        style = "solid"
        if u == v:  # self-loop = recurrent link
            style = "dotted"
        ax annotate("", xy(pos[v], xycoordinates="data"),
                   xytext(pos[u], textcoordinates="data"),
                   arrowprops=dict(arrowstyle="->",
                                   linestyle=style,
                                   linewidth=width,
                                   color=color,
                                   alpha=alpha,
                                   shrinkA=10, shrinkB=10)))

    plt.title(title)
```
92
93
94
95
96
97
98
99
100
101
102
103

`plt.axis("off")`
`plt.tight_layout()`
`if save_as:`
    `plt.savefig save_as, dπ=300, bbox_inches="tight")`
`plt.show()`

`# ----------------------------`
`# Example run`
`# ----------------------------`

`if __name__ == "__main__":`
    `G = generate_gnarl_genome()`
    `plot_gnarl(G, title="GNARL Example Genome", save_as="gnarl_genome.png")`

### 16.5.5 Neuroevolution of augmenting topologies (NEAT)

The Neuroevolution of Augmenting Topologies (NEAT) method is a powerful approach in artificial intelligence that applies evolutionary algorithms to optimize both the weights and the topology of artificial neural networks (ANNs). Traditional neuroevolution methods generally evolve fixed-topology networks, limiting their ability to discover innovative structures. NEAT circumvents this limitation by dynamically augmenting network topology, thereby evolving increasingly complex structures while maintaining genetic compatibility through historical markings. Let the genome representation of a neural network be denoted as a directed graph $G = (V, E)$, where the set $V$ consists of input, hidden, and output neurons, and $E$ represents weighted connections between these neurons. Each connection has an associated weight $w_{ij}$ such that the activation $a_j$ of neuron $j$ is given by

$$a_j = f \left( \sum_{i \in \mathcal{P}(j)} w_{ij} a_i + b_j \right) \qquad (16.394)$$

where $\mathcal{P}(j)$ is the set of neurons feeding into neuron $j$, $b_j$ is the bias term, and $f(\cdot)$ is the activation function, which is often chosen to be a non-linear function such as the sigmoid

$$f(x) = \frac{1}{1 + e^{-x}} \qquad (16.395)$$

or a rectified linear unit (ReLU)

$$f(x) = \max(0, x). \qquad (16.396)$$

In NEAT, the evolutionary process begins with a population of simple neural networks (often with only input and output layers) and progressively adds new connections and new neurons via mutation operations. The weight evolution follows a traditional genetic algorithm, where mutations modify the connection strengths:

$$w_{ij}^{(t+1)} = w_{ij}^{(t)} + \eta \cdot \mathcal{N}(0, \sigma^2) \qquad (16.397)$$

where $\eta$ is the learning rate and $\mathcal{N}(0, \sigma^2)$ is a Gaussian noise term. Structural mutations can add a new connection between two previously unconnected nodes with a probability $p_c$, or a new node can be inserted into an existing connection with probability $p_n$. The insertion of a new node $v_k$ replaces an edge $(i, j)$ with two edges $(i, k)$ and $(k, j)$, and their weights are initialized as

$$w_{ik}^{(0)} = w_{ij}, \quad w_{kj}^{(0)} = 1. \qquad (16.398)$$

One of the central innovations of NEAT is the use of historical markings to track the lineage of genes (connections). Each connection gene is assigned an innovation number upon creation, which
<figure><img src="image_217.png" /><figcaption>Figure 16.48: Visualization of a NEAT genome: Input nodes (blue), hidden nodes (green), and output nodes (red) represent the network structure. Connections are colored blue for positive weights and red for negative weights, with thickness proportional to weight magnitude. Enabled connections are shown as solid arrows, while disabled ones are dashed</figcaption></figure>

remains unchanged throughout evolution. The similarity $d$ between two genomes $G_1$ and $G_2$ is computed using the genetic distance metric

$$d = \frac{c_1 E}{N} + \frac{c_2 D}{N} + c_3 W \qquad (16.399)$$

where $E$ is the number of excess genes, $D$ is the number of disjoint genes, $W$ is the average weight difference of matching genes, and $c_1, c_2, c_3$ are scaling coefficients. The normalization factor $N$ accounts for genome length to prevent excessive penalties for larger networks. This metric is crucial for implementing speciation, where similar networks are grouped into species to encourage innovation without being prematurely eliminated by competition. The fitness of each network $G$ is evaluated using a task-dependent function $F(G)$, and species-level fitness sharing is applied to maintain diversity:

$$F'_{i} = \frac{F_{i}}{\sum_{j \in \text{species}} s(d_{ij})} \qquad (16.400)$$

where $s(d_{ij})$ is a step function that determines whether two genomes belong to the same species:

$$s(d_{ij}) = \begin{cases} 1, & \text{if } d_{ij} < \delta \\ 0, & \text{otherwise} \end{cases} \qquad (16.401)$$

with $\delta$ being a speciation threshold. Networks within a species compete primarily amongst themselves, fostering the survival of innovative structures. The crossover operation between two parent
genomes $G_1$ and $G_2$ is performed by inheriting matching genes randomly and retaining excess and disjoint genes from the more fit parent:

$$G_{\text{child}} = \{g_i | g_i \in G_1 \cap G_2 \text{ with probability } 0.5\} \cup \{g_i | g_i \in G_{\text{fitter}}\} \quad (16.402)$$

where $G_{\text{fitter}}$ is the genome with the higher fitness. If a connection gene in one parent is disabled due to mutation, it is inherited in a disabled state unless reactivated by future mutations. Another fundamental aspect of NEAT is complexification, wherein networks start minimally and gradually increase in complexity by adding nodes and connections. Unlike fixed-topology methods, NEAT allows evolution to discover increasingly sophisticated representations, leading to efficient problem-solving. The network complexity at generation $t$ is measured by

$$C(t) = |V| + |E| \quad (16.403)$$

where $|V|$ and $|E|$ denote the number of neurons and connections, respectively. Over time, the structural complexity increases as new mutations introduce novel architectures. In practice, the evolutionary process iterates until convergence, defined by an optimal fitness threshold $F_{\text{opt}}$ such that

$$\max_{G \in \text{population}} F(G) \geq F_{\text{opt}}. \quad (16.404)$$

The mutation rates $p_c$ and $p_n$ are typically annealed to balance exploration and exploitation, governed by

$$p_c(t) = p_c(0)e^{-\lambda t}, \quad p_n(t) = p_n(0)e^{-\mu t} \quad (16.405)$$

where $\lambda, \mu$ are decay rates. In summary, NEAT dynamically evolves both weights and network structures while preserving historical innovations and maintaining diversity through speciation. The key principles—incremental complexification, historical markings, and speciation—enable NEAT to efficiently discover topologies that traditional methods struggle to optimize. The continuous balance between exploration and exploitation, governed by well-defined evolutionary operators, results in a robust and self-adaptive neuroevolution framework.

### 16.5.5.1 Python Code to generate Figure 16.48

The NeuroEvolution of Augmenting Topologies (NEAT) algorithm evolves both the weights and topology of neural networks. A plot for NEAT should show:

1. Input, hidden, and output nodes.

2. Edges with weights (positive = blue, negative = red).

3. Disabled connections (often drawn as dashed or faded).

4. Labels showing node IDs.

The Python code below produces the Figure 16.48 illustrating Neuroevolution of augmenting topologies (NEAT).

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

```python
import networkx as nx
```

```python
import random
```

```python
# -----------------------------
```

```python
# Example NEAT genome
```

```python
# -----------------------------
```

```python
def generate_neat_genome(numInputs=3, numHidden=4, numOutputs=2, seed
    =42):
```

```python
    random.seed(seed)
```
```python
G = nx.DiGraph()

# Input nodes (layer = 0)
inputs = [f"I{i}" for i in range(1, numInputs+1)]
for i, n in enumerate inputs:
    G.add_node(n, layer="input", pos=(i, 0))

# Hidden nodes (layer = 1)
hidden = [f"H{i}" for i in range(1, numHidden+1)]
for i, n in enumerate(hidden):
    x, y = random.random(0, numInputs-1), random.random(0.5, 1.5)
    G.add_node(n, layer="hidden", pos=(x, y))

# Output nodes (layer = 2)
outputs = [f"O{i}" for i in range(1, numOutputs+1)]
for i, n in enumerate outputs:
    G.add_node(n, layer="output", pos=(i, 2))

# Add edges with weights
edges = []
for inp in inputs:
    for hid in hidden:
        if random.random() < 0.6:
            edges.append((inp, hid))
    for out in outputs:
        if random.random() < 0.7:
            edges.append((hid, out))
    for inp in inputs:
        for out in outputs:
            if random.random() < 0.3:
                edges.append((inp, out))

for (u, v) in edges:
    G.add_edge(u, v,
               weight=round(random.random()(-1, 1), 2),
               enabled=(random.random() > 0.2)  # some connections
                                   disabled

return G

# ------------------------
# Plotting function
# ------------------------
def plot_neat(G, title="NEAT Genome", save_as=None):
    pos = {n: d["pos"] for n, d in G.nodes(data=True)}

# Node colors
color_map = []
for n, d in G.nodes(data=True):
    if d["layer"] == "input":
        color_map.append("lightblue")
    elif d["layer"] == "hidden":
        color_map.append("lightgreen")
    else:
        return G
```
```python
color_map.append("salmon")

# Draw nodes
nx.draw_networkx_nodes(G, pos, node_color color_map, node_size=800,
                       edgecolors="k")
nx.draw_networkx labels(G, pos, font_size=8)

# Draw edges
ax = plt.gca()
for u, v, d in GEdges(data=True):
    w = d["weight"]
    enabled = d["enabled"]
    color = "blue" if w > 0 else "red"
    width = 1 + 2 * abs(w)
    style = "solid" if enabled else "dashed"

    ax annotate("", xy(pos[v], xycoordinates="data"),
              xytext(pos[u], textcoordinates="data"),
              arrowprops=dict(arrowstyle="->",
                              style=linestyle,
                              linewidth=width,
                              color=color,
                              alpha=0.6,
                              shrinkA=10, shrinkB=10))

    plt.title(title)
    plt.axis("off")
    plt.tight_layout()
    if save_as:
        plt.savefig save_as, dπ=300, bbox_inches="tight")
    plt.show()

# Example run
if __name__ == "__main__":
    G = generate_neat_genome()
    plot_neat(G, title="NEAT Example Genome", save_as="neat_genome.png")
```

### 16.5.6 Hypercube-based NeuroEvolution of Augmenting Topologies (HyperNEAT)

The Hypercube-based NeuroEvolution of Augmenting Topologies (HyperNEAT) method is a neu-
roevolutionary approach that extends the NeuroEvolution of Augmenting Topologies (NEAT) al-
gorithm by leveraging a Compositional Pattern Producing Network (CPPN) to encode the con-
nectivity patterns of artificial neural networks (ANNs) through a substrate. The fundamental
idea behind HyperNEAT is to exploit geometric regularities in the problem domain by encoding
the connectivity of a neural network as a function of spatial coordinates, allowing it to generate
large-scale neural structures with intricate connectivity patterns in a computationally efficient and
evolutionarily scalable manner. Mathematically, let the neural substrate be defined as a set of nodes
embedded in a D-dimensional Euclidean space, where each node is assigned a spatial coordinate
x ∈ R<sup>D</sup>. The weight matrix W of the neural network is not directly evolved but rather generated
<figure><img src="image_218.png" /><figcaption>Figure 16.49: Visualization of a HyperNEAT substrate-based neural network: Input nodes (blue), hidden nodes (green), and output nodes (red) are positioned in a geometric substrate. Edge thickness is proportional to weight magnitude, with blue representing positive weights and red representing negative weights</figcaption></figure>

dynamically using a function $f: \mathbb{R}^{2D} \to \mathbb{R}$ that maps the spatial positions of node pairs to connection weights. This function $f$ is instantiated by a Compositional Pattern Producing Network (CPPN), which is a fully connected feedforward neural network with activation functions that can include sigmoidal, Gaussian, sine, and identity functions, among others. The connectivity weight between two neurons located at coordinates $\mathbf{x}_i$ and $\mathbf{x}_j$ is given by:

$$W_{ij} = f(\mathbf{x}_i, \mathbf{x}_j) \qquad (16.406)$$

where the function $f(\mathbf{x}_i, \mathbf{x}_j)$ is parameterized by the **weights and activation functions** of the CPPN, which itself undergoes **neuroevolution** using the NEAT algorithm. The CPPN represents a **compressed encoding** of the connectivity pattern, which allows the network to **exploit spatial regularities** in a way that would be infeasible with a direct encoding approach. The NEAT algorithm evolves the CPPN by **mutating and recombining genomes**, where each genome $G$ represents a **set of neurons** and **connections** forming the CPPN. A genome $G$ in the NEAT representation consists of a **set of node genes** and a **set of connection genes**, where each connection gene is represented as a tuple:

$$C_{ij} = (i, j, w_{ij}, \text{enabled}) \qquad (16.407)$$

where $i$ and $j$ are node indices, $w_{ij}$ is the connection weight, and enabled is a binary flag indicating whether the connection is active. The fitness function $F(G)$ that guides the evolution is determined by the performance of the **decoded ANN** on the given task. The evolution proceeds using **mutation operators** such as:

$$w_{ij} \leftarrow w_{ij} + \mathcal{N}(0, \sigma^2) \qquad (16.408)$$
where $\mathcal{N}(0, \sigma^2)$ is a **Gaussian perturbation**, and structural mutations that **add nodes and connections**, preserving the historical marking mechanism of NEAT. Once the CPPN is evolved, it is queried at every pair of spatial coordinates $(\mathbf{x}_i, \mathbf{x}_j)$ in the neural substrate to **generate the weight matrix**:

$$W = \{f(\mathbf{x}_i, \mathbf{x}_j) \mid \forall (\mathbf{x}_i, \mathbf{x}_j) \in S \times S\} \qquad (16.409)$$

where $S$ is the set of nodes in the substrate. This procedure ensures that the evolved ANN **inherits the geometric properties** encoded by the CPPN, leading to **topological coherence** and the ability to **scale up** efficiently. The substrate structure can be fixed or evolved, leading to a variety of architectures including **multi-layer perceptrons**, **convolutional-like structures**, and **even more complex topologies**. The implicit encoding provided by the CPPN enables **patterned connectivity** with **symmetries**, **repetitions**, and **other motifs** that are **biologically plausible** and **computationally advantageous**. HyperNEAT's search space is thus fundamentally different from traditional neuroevolution approaches, as it searches over **functions that encode networks** rather than over networks themselves. This means that small changes in the CPPN can lead to **large, structured modifications** in the ANN's connectivity, a property known as **indirect encoding**. Formally, the space of networks $N$ encoded by a CPPN is given by:

$$N = \{f(\mathbf{x}_i, \mathbf{x}_j) \mid \mathbf{x}_i, \mathbf{x}_j \in \mathbb{R}^D\} \qquad (16.410)$$

where $f$ is constrained by the activation functions and weight structure of the CPPN. The expressiveness of the encoding is controlled by the **activation functions** in the CPPN, which determine the nature of the patterns it can represent. A significant advantage of HyperNEAT is its ability to **exploit domain regularities** by leveraging **geometric relationships** within the substrate. For example, in a **vision-based task**, neurons representing nearby pixels should have stronger connections, which can be captured by a **distance-sensitive** function:

$$f(\mathbf{x}_i, \mathbf{x}_j) = \alpha \exp(-\beta \|\mathbf{x}_i - \mathbf{x}_j\|^2) \qquad (16.411)$$

where $\alpha$ and $\beta$ are **evolved parameters** that control the strength and scale of connectivity. This equation ensures that neurons **form meaningful topologies** based on spatial structure. HyperNEAT optimizes the **encoding function** rather than individual weights, allowing it to **generalize connectivity patterns** across different domains and tasks. This ability to evolve **topological regularities** has been successfully applied in **robot control**, **game playing**, **function approximation**, and **large-scale pattern recognition**.

### 16.5.6.1 Python Code to generate Figure 16.49

HyperNEAT generates neural networks by mapping a geometric substrate via a CPPN.

1. Nodes are in a 2D (or 3D) substrate.

2. Connections are determined by the CPPN, usually sparse but patterned.

3. Visualization should highlight:

(a) Node positions (substrate coordinates)

(b) Connection weights (width proportional, color by sign)

The Python code below produces the Figure 16.49 illustrating Hypercube-based NeuroEvolution of Augmenting Topologies (HyperNEAT).

1. import matplotlib.pyplot as plt

2. import networkx as nx

3. import random

4.
```python
5 # -----------------------------
6 # Example HyperNEAT substrate
7 # -----------------------------
8 def generate_hyperneat_substrate(numInputs=4, num_hidden=6, numOutputs
9                                     =2, seed=42):
10         random.seed(seed)
11         G = nx.DiGraph()
12
13         # Input nodes (y=0)
14         inputs = [f"I{i}" for i in range(1, numInputs+1)]
15         for i, n in enumerate(Inputs):
16             G.add_node(n, layer='input', pos=(i,0), activation='linear')
17
18         # Hidden nodes (y(random)
19         hidden = [f"H{i}" for i in range(1, num_hidden+1)]
20         for i, n in enumerate(hidden):
21             x, y = randomhea
22             act = random.choice([ 'tanh', 'sigmoid', 'relu'])
23             G.add_node(n, layer='hidden', pos=(x,y), activation=act)
24
25         # Output nodes (y=2)
26         outputs = [f"O{i}" for i in range(1, numOutputs+1)]
27         for i, n in enumerate(Outputs):
28             G.add_node(n, layer='output', pos=(i,2), activation='tanh')
29
30         # Connect input to hidden
31         for inp in inputs:
32             for hid in hidden:
33                 if random.random() < 0.7:
34                     G.add_edge inp, hid, weight=round(randomhea
35                 )
36
37         # Connect hidden to output
38         for hid in hidden:
39             for out in outputs:
40                 if random.random() < 0.8:
41                     G.add_edge hid, out, weight=round(randomhea
42                 )
43         return G
44
# -----------------------------
45 # Plotting function
46 # -----------------------------
47 def plot_hyperneat(G, title="HyperNEAT Substrate", save_as=None):
48         pos = {n[ 'pos'] for n, data in G.nodes(data=True)}
49
50         # Node colors
51         color_map = []
52         for n, d in G.nodes(data=True):
53             if d['layer'] == 'input': color_map.append('lightblue')
54             elif d['layer'] == 'hidden': color_map.append('lightgreen')
55             else: color_map.append('salmon')
56
57         # Draw nodes
58         nx.draw_networkx_nodes(G, pos, node_color=co
59         edgecolors='k')
```
```
nx.draw_networkx labels(G, pos, labels={n: f"{n}\n{d['activation']}"}
    for n, d in G.nodes(data=True)}, font_size=8)

# Draw edges
ax = plt.gca()
for u, v, d in GEdges(data=True):
    w = d['weight']
    color = 'blue' if w>0 else 'red'
    width = 1 + 2*abs(w)
    ax annotate('',
        xy(pos[v], xycoordinates='data',
            xytext(pos[u], textcoordinates='data',
                arrowprops(dict(arrowstyle='->', linewidth=width,
                    color=color, alpha=0.6,
                        shrinkA=10, shrinkB=10))
            ))

    plt.title(title)
    plt.axis('off')
    plt.tight_layout()
    if save_as: plt.savefig save_as, dπ=300, bbox_inches='tight')
    plt.show()

# Example run
if __name__ =="__main__":
    G = generate_hyperneat_substrate()
    plot_hyperneat(G, title="HyperNEAT Example Genome", save_as="
        hyperneat_genome.png")
```

### 16.5.7 Evolvable Substrate Hypercube-based NeuroEvolution of Augmenting Topologies (ES-HyperNEAT)

The Evolvable Substrate Hypercube-based NeuroEvolution of Augmenting Topologies (ES-HyperNEAT) method is an advanced extension of the Hypercube-based NeuroEvolution of Augmenting Topologies (HyperNEAT) algorithm, which itself is based on NeuroEvolution of Augmenting Topologies (NEAT). The fundamental premise of ES-HyperNEAT is the encoding of connectivity patterns in neural networks using Compositional Pattern-Producing Networks (CPPNs), which are evolved through genetic algorithms. This approach allows for the indirect encoding of large-scale artificial neural networks (ANNs) by leveraging geometric regularities and adaptive connectivity growth to discover task-relevant neural substrates in a computationally efficient manner. The ES-HyperNEAT methodology extends HyperNEAT by **evolving connectivity patterns based on local fitness evaluations within a neural substrate**, rather than relying on global constraints alone. The key innovation in ES-HyperNEAT is the ability to discover and refine important network regions dynamically, without explicitly requiring a predefined topology. Given a **substrate space $S$** with spatial coordinates for nodes $x \in \mathbb{R}^d$, the connection strength between two nodes $x_i, x_j \in S$ is encoded by a **functionally evolved CPPN**, represented mathematically as

$$w_{ij} = f_{\text{CPPN}}(\mathbf{x}_i, \mathbf{x}_j). \qquad (16.412)$$

Here, $f_{\text{CPPN}}$ is a neural network that encodes the weights between nodes based on their spatial location. The CPPN is evolved via **NEAT**, which optimizes its parameters and topology to generate increasingly complex connectivity patterns. Unlike conventional ANN weight optimization,
<figure><img src="image_219.png" /><figcaption>Figure 16.50: Visualization of an ES-HyperNEAT substrate-based neural network: Input nodes (blue), hidden nodes (green), and output nodes (red) are positioned in a geometric substrate. Edge thickness is proportional to weight magnitude, with blue representing positive weights and red representing negative weights</figcaption></figure>

which directly evolves individual weights, ES-HyperNEAT leverages the implicit regularities in the substrate to determine connections based on geometric principles. To formally describe the **weight expression process**, let us define the CPPN-generated weight function as a mapping:

$$W : \mathbb{S} \times \mathbb{S} \to \mathbb{R}, \quad (\mathbf{x}_i, \mathbf{x}_j) \mapsto w_{ij}. \qquad (16.413)$$

A fundamental property of this mapping is **substrate symmetry**, which is captured by the fact that weight magnitudes for symmetric inputs in the domain follow similar functional properties:

$$\forall (\mathbf{x}_i, \mathbf{x}_j) \in \mathbb{S} \times \mathbb{S}, \quad W(\mathbf{x}_i, \mathbf{x}_j) \approx W(\mathbf{x}_j, \mathbf{x}_i). \qquad (16.414)$$

This symmetry allows for computational efficiency, as only a fraction of the connections need to be explicitly evaluated. In ES-HyperNEAT, unlike conventional HyperNEAT, the adaptive growth mechanism is introduced by defining an adaptive threshold function $\tau(\mathbf{x})$ over the substrate, which dynamically determines whether a connection should be expressed. The connection viability is determined by:

$$|w_{ij}| > \tau(\mathbf{x}_i, \mathbf{x}_j), \qquad (16.415)$$

where $\tau$ is a learned function that varies over the substrate space and is itself evolved during the learning process. The introduction of adaptive thresholds enables ES-HyperNEAT to selectively prune redundant connections while preserving the essential connectivity for task-relevant computation. The neural network instantiated from the CPPN-generated substrate follows a typical activation dynamics governed by the weighted sum and activation function:

$$a_i = \sigma \left( \sum_j w_{ij} a_j + b_i \right), \qquad (16.416)$$
where $a_i$ is the activation of node $i$, $w_{ij}$ is the synaptic weight from node $j$ to node $i$, $b_i$ is a bias term, and $\sigma$ is a nonlinear activation function (e.g., sigmoid, tanh, or ReLU). The ES-HyperNEAT approach further refines the weight refinement strategy by incorporating a local function-based weight refinement, denoted as $W'$, which modifies the initial weights $W$ using a learned **adaptive function** $g(\mathbf{x}_i, \mathbf{x}_j)$:

$$w'_{ij} = W(\mathbf{x}_i, \mathbf{x}_j) + g(\mathbf{x}_i, \mathbf{x}_j). \qquad (16.417)$$

This refinement process iteratively enhances connectivity by **reinforcing important regions** in the substrate while **eliminating weak connections** through thresholding and weight decay. The effect is a progressively **optimized connectivity pattern** that dynamically adapts to the learning environment. From a mathematical perspective, the evolution of the CPPN itself follows **genetic mutation and crossover dynamics**, where the genetic encoding of the CPPN, denoted by $\Theta_{CPPN}$, undergoes variation through mutation functions $\mathcal{M}$ and crossover functions $\mathcal{C}$:

$$\Theta_{CPPN}^{(t+1)} = \mathcal{M}(\mathcal{C}(\Theta_{CPPN}^{(t)}, \Theta_{CPPN}^{(t')})). \qquad (16.418)$$

Here, $t$ represents the generation index, and $\Theta_{CPPN}^{(t)}$ and $\Theta_{CPPN}^{(t')}$ are the CPPN parameters of parent solutions selected based on fitness evaluations. The fitness function $F$ evaluates the network's performance on a given task and is defined as:

$$F(\Theta_{CPPN}) = \sum_{k} \mathcal{L}(\hat{y}_k, y_k), \qquad (16.419)$$

where $\mathcal{L}$ is a loss function comparing predicted outputs $\hat{y}_k$ to true labels $y_k$. The evolutionary process seeks to maximize $F$, leading to increasingly optimized neural architectures. Overall, the **ES-HyperNEAT** approach represents a mathematically rigorous framework for evolving complex neural network topologies through **adaptive connectivity growth**, **thresholded weight refinement**, and **neurogenesis-driven substrate expansion**. The resulting architectures are not only **topologically optimized** but also **computationally scalable**, making ES-HyperNEAT a powerful tool for the evolution of high-dimensional neural representations in artificial intelligence.

### 16.5.7.1 Python Code to generate Figure 16.50

ES-HyperNEAT generates large-scale neural networks by evolving substrates and connection patterns via Compositional Pattern Producing Networks (CPPNs).

1. Nodes are arranged in a geometric substrate (usually 2D or 3D).

2. Connections are created based on the CPPN output across node coordinates.

3. Visualizations typically show:

4. Node positions in substrate $(x, y)$.

5. Weighted connections (colored by sign, thickness by magnitude).

The Python code below produces the Figure 16.51 illustrating Evolvable Substrate Hypercube-based NeuroEvolution of Augmenting Topologies (ES-HyperNEAT).

```python
import matplotlib.pyplot as plt
import networkx as nx
import random

# Example ES-HyperNEAT substrate
```
```python
def generate_eshyperneat_substrate(numInputs=4, numHidden=6,
                                     numOutputs=2, seed=42):
    random.seed(seed)
    G = nx.DiGraph()

    # Input nodes (positioned at y=0)
    inputs = [f"I{i}" for i in range(1, numInputs+1)]
    for i, n in enumerate(Inputs):
        G.add_node(n, layer='input', pos=(i, 0), activation='linear')

    # Hidden nodes (random positions in middle y-layer)
    hidden = [f"H{i}" for i in range(1, numHidden+1)]
    for i, n in enumerate(hidden):
        x, y = randomhea
        act = random.choice([ 'tanh', 'sigmoid', 'relu'])
        G.add_node(n, layer='hidden', pos=(x, y), activation=act)

    # Output nodes (y=2)
    outputs = [f"O{i}" for i in range(1, numOutputs+1)]
    for i, n in enumerate(Outputs):
        G.add_node(n, layer='output', pos=(i, 2), activation='tanh')

    # Connect input to hidden
    for inp in Inputs:
        for hid in Hidden:
            if random.random() < 0.7:
                G.add_edge inp, hid, weight=round(randomhea
                                                (-1, 1), 2))

    # Connect hidden to output
    for hid in Hidden:
        for out in Outputs:
            if random.random() < 0.8:
                G.add_edge(hid, out, weight=round(randomhea
                                                (-1, 1), 2))

    return G
```

```python
# -----------------------------
```

```python
# Plotting function
```

```python
# -----------------------------
```

```python
def plot_eshyperneat(G, title="ES-HyperNEAT Substrate", save_as=None):
    pos = {n:sdata[ 'pos'] for n, data in G.nodes(data=True)}
```

```python
# Node colors
```

```python
color_map = []
```

```python
for n, d in G.nodes(data=True):
    if d['layer'] == 'input': color_map.append('lightblue')
    elif d['layer'] == 'hidden': color_map.append('lightgreen')
    else: color_map.append('salmon')
```

```python
# Draw nodes
```

```python
nx.draw_networkx_nodes(G, pos, node_color=co
edgecolors='k')
```

```python
nx.draw_networkx labels(G, pos, labels={n:f"{n} 
"d['activation']}")
```
```python
# Draw edges
ax = plt.gca()
for u,v,d in GEdges(data=True):
    w = d['weight']
    color = 'blue' if w>0 else 'red'
    width = 1 + 2*abs(w)
    ax annotate('',
        xy(pos[v], xycoordinates='data',
        xytext(pos[u], textcoordinates='data',
        arrowprops(dict(arrowstyle='->', linewidth=width,
            color=color, alpha=0.6,
            shrinkA=10, shrinkB=10))
        ),
        title(title)
        axis('off')
        tight_layout()
        if save_as: plt.savefig save_as, dπ=300, bbox_inches='tight')
        show()

# Example run
if __name__ == "__main__":
    G = generate_eshyperneat_substrate()
    plot_eshyperneat(G, title="ES-HyperNEAT Example Genome", save_as="eshyperneat_genome.png")
```

## 16.5.8 Evolutionary Acquisition of Neural Topologies (EANT/EANT2)

**Literature Review:** EANT has been assessed on multiple benchmark tasks, such as the RoboCup keepaway challenge (Metzen et. al. (2007) [986]) and the double-pole balancing problem (Kassahun and Sommer 2005 [987]), consistently achieving high performance. A more advanced version, EANT2, was later introduced and tested on a visual servoing task, where it surpassed both the traditional iterative Gauss-Newton method and NEAT (Siebel and Sommer [988]). Additional experiments have also been conducted to evaluate its capabilities in classification problems (Siebel and Sommer [990]).

**Recent Literature Review:**

The Evolutionary Acquisition of Neural Topologies (EANT) and its extension EANT2 are advanced methods for evolving artificial neural networks (ANNs) through an evolutionary process that simultaneously optimizes both the topology and the weights of the network. The primary principle underlying EANT/EANT2 is the combination of evolutionary algorithms (EAs) and neural network training, allowing for the automatic discovery of optimal network architectures without the need for manual design. The method begins with an initial minimal structure and progressively complexifies through mutations and recombinations, ensuring efficient exploration of the search space while avoiding unnecessary complexity. A neural network in the context of EANT/EANT2 can be defined as a directed graph $G = (V, E)$, where the set of vertices $V$ represents neurons, and the set of edges $E$ represents synaptic connections. Each edge $e_{ij} \in E$ carries an associated weight $w_{ij}$, which modulates the signal transmission from neuron $i$ to neuron $j$. Mathematically,
the activation $a_j$ of neuron $j$ at time step $t$ is given by:

$$a_j(t) = f\left(\sum_{i \in \text{Pre}(j)} w_{ij} a_i(t-1) + b_j\right), \qquad (16.420)$$

where Pre($j$) denotes the set of presynaptic neurons to $j$, $w_{ij}$ are the connection weights, $b_j$ is the bias term, and $f(\cdot)$ is the activation function, typically chosen as the sigmoid function

$$f(x) = \frac{1}{1 + e^{-\lambda x}} \qquad (16.421)$$

or the hyperbolic tangent function

$$f(x) = \tanh(x). \qquad (16.422)$$

In the EANT framework, network evolution begins with a set of simple neural networks containing only input and output neurons, with hidden neurons introduced incrementally via evolutionary operators. The evolution of network topology is governed by genetic operators such as mutation

<figure><img src="image_220.png" /><figcaption>Figure 16.51: Visualization of the Evolutionary Acquisition of Neural Topologies (EANT2): Input nodes (blue), hidden nodes (green), and output nodes (red) form an evolved topology. Edge colors represent weights (blue: positive, red: negative) with thickness proportional to magnitude</figcaption></figure>

and crossover. Mutation can involve perturbing connection weights, adding or deleting neurons, or modifying the topology in a more structured manner. The mutation of a weight follows a Gaussian perturbation model:

$$w_{ij}^{(t+1)} = w_{ij}^{(t)} + \eta \mathcal{N}(0, \sigma^2), \qquad (16.423)$$

where $\mathcal{N}(0, \sigma^2)$ is a Gaussian-distributed random variable with mean zero and variance $\sigma^2$, and $\eta$ is a learning rate parameter. The selection process follows a fitness-based evaluation, where each neural network $N_i$ in a population $\mathcal{P} = \{N_1, N_2, \dots, N_M\}$ is assigned a fitness score $F(N_i)$ based
on its performance on a given task. The selection probability of an individual network follows a Boltzmann distribution:

$$P(N_i) = \frac{e^{\beta F(N_i)}}{\sum_{j=1}^{M} e^{\beta F(N_j)}}, \qquad (16.424)$$

where $\beta$ is a temperature parameter controlling selection pressure. EANT2 extends EANT by introducing covariance matrix adaptation (CMA) for more efficient weight optimization. In CMA, the evolution of the weight vector $\mathbf{w}$ is guided by a covariance matrix $\mathbf{C}$, which dynamically adapts to the topology of the search landscape. The weight update rule in EANT2 is given by:

$$\mathbf{w}^{(t+1)} = \mathbf{w}^{(t)} + \mathbf{C}^{1/2}\mathcal{N}(\mathbf{0}, \mathbf{I}), \qquad (16.425)$$

where $\mathbf{I}$ is the identity matrix. The covariance matrix is updated using an evolution path formulation:

$$\mathbf{C}^{(t+1)} = (1 - c_c)\mathbf{C}^{(t)} + c_c\mathbf{p}_c\mathbf{p}_c^\top, \qquad (16.426)$$

where $c_c$ is a learning rate, and $\mathbf{p}_c$ is the evolution path. Throughout the evolutionary process, EANT/EANT2 ensures that only structurally beneficial changes are retained, leading to an efficient exploration-exploitation tradeoff. The incremental growth of topology in EANT can be mathematically formalized by defining the probability of adding a new neuron $v_k$ as

$$P_{\text{add}}(v_k) = \frac{1}{1 + e^{-\gamma \Delta F}}, \qquad (16.427)$$

where $\Delta F$ represents the fitness improvement from the structural modification, and $\gamma$ is a control parameter. At each iteration, the population of neural networks undergoes selection, mutation, and recombination, leading to the formation of an improved generation. The optimization objective is to maximize a performance metric $J(\mathbf{w}, \mathcal{T})$, where $\mathcal{T}$ represents the training dataset. The overall evolutionary update rule can be expressed as

$$\mathbf{w}^{(t+1)} = \mathbf{w}^{(t)} + \alpha \nabla_{\mathbf{w}} J(\mathbf{w}, \mathcal{T}) + \mathbf{C}^{1/2} \mathcal{N}(\mathbf{0}, \mathbf{I}), \qquad (16.428)$$

where $\alpha$ is a learning rate and the second term represents stochastic exploration. The convergence of EANT/EANT2 is governed by the stability of the weight adaptation dynamics. If the eigenvalues of the Jacobian matrix

$$\mathbf{J} = \frac{\partial f(\mathbf{w})}{\partial \mathbf{w}} \qquad (16.429)$$

lie within the unit circle, the evolutionary process stabilizes, ensuring convergence. Otherwise, further structural modifications are required to regularize the topology. Thus, EANT/EANT2 provides a mathematically rigorous framework for evolving neural networks by simultaneously optimizing topology and weights, leveraging evolutionary principles, and incorporating efficient weight adaptation techniques such as CMA. The fundamental strength of the method lies in its ability to construct minimal yet powerful architectures that efficiently learn complex functions.

### 16.5.8.1 Python Code to generate Figure 16.51

The Evolutionary Acquisition of Neural Topologies (EANT) algorithm evolves neural network topologies using evolutionary strategies. A plot should show:

1. Input layer, hidden nodes, output layer.

2. Evolved connections (weighted, possibly sparse).

3. Mutations such as added nodes or connections.

4. Generation info (to highlight evolutionary growth).
The Python code below produces the Figure 16.51 illustrating Evolutionary Acquisition of Neural Topologies (EANT) algorithm.

```python
import networkx as nx
import matplotlib.pyplot as plt
import random

def generate_eant_genome(numInputs=3, num_hidden=5, numOutputs=2, seed
=42):
    random.seed(seed)
    G = nx.DiGraph()

    # Input layer
    inputs = [f"I{i}" for i in range(1, numInputs + 1)]
    for n in inputs:
        G.add_node(n, layer="input", activation="linear")

    # Hidden layer (evolvable)
    hidden = [f"H{i}" for i in range(1, num_hidden + 1)]
    for n in hidden:
        act = random.choice(["tanh", "relu", "sigmoid"])
        G.add_node(n, layer="hidden", activation=act)

    # Output layer
    outputs = [f"O{i}" for i in range(1, numOutputs + 1)]
    for n in outputs:
        G.add_node(n, layer="output", activation="softmax")

    # Connect inputs to hidden
    for i in inputs:
        for h in hidden:
            if random.random() < 0.6:
                G.add_edge(i, h, weight=round(randomhea
n(-1, 1), 2))

    # Connect hidden to hidden (possible recurrent/self connections)
    for h1 in hidden:
        for h2 in hidden:
            if h1 != h2 and random.random() < 0.2:
                G.add_edge(h1, h2, weight=round(random
n(-1, 1), 2))

    # Connect hidden to outputs
    for h in hidden:
        for o in outputs:
            if random.random() < 0.7:
                G.add_edge(h, o, weight=round(random
n(-1, 1), 2))

    return G

def plot_eant(G, title="EANT Genome Example", save_as=None):
    # Node colors
    color_map = []
    for n, data in G.nodes(data=True):
        if data["layer"] == "input":
```
```python
color_map.append("lightblue")
elif data["layer"] == "hidden":
    color_map.append("lightgreen")
else:
    color_map.append("salmon")

# Positions by layer
pos = {}
inputs = [n for n, d in G.nodes(data=True) if d["layer"] == "input"]
hidden = [n for n, d in G.nodes(data=True) if d["layer"] == "hidden"]
outputs = [n for n, d in G.nodes(data=True) if d["layer"] == "output"]

for i, n in enumerate(Inputs):
    pos[n] = (-1, i)

for i, n in enumerate(hidden):
    pos[n] = (0, i)

for i, n in enumerate(Outputs):
    pos[n] = (1, i)

# Draw graph
plt.figure(figsize=(10, 6))
nx.draw_networkx_nodes(G, pos, node_color=Color_map, node_size=800,
                       edgecolors="black")
nx.draw_networkx labels(G, pos, font_size=8)

# Weighted edges
for u, v, d in GEdges(data=True):
    w = d["weight"]
    color = "blue" if w > 0 else "red"
    width = 1 + abs(w) * 2
    alpha = 0.6
    nx.draw_networkx_edges(G, pos, edgelist=[(u, v)], width=width,
                           edge_color=Color, alpha=alpha)

# Show weights
edge_labels = { (u, v): d["weight"] for u, v, d in G.Edges(data=True) }
nx.draw_networkx edge labels(G, pos, edge labels=edge labels,
                             font size=7)

plt. title (title)
plt. axis ("off")
if save_as:
    plt. savefig (save_as, dpi=300, bbox_inches="tight")
    plt. show()

if __name__ == "__main__":
    G = generate_eant_genome()
    plot_eant(G, title="Evolutionary Acquisition of Neural Topologies (
                       EANT2)", save_as="eant2_genome.png")
```
### 16.5.9 Interactively Constrained Neuro-Evolution (ICONE)

The Interactively Constrained Neuro-Evolution (ICONE) method in artificial intelligence is a sophisticated approach that integrates interactive constraints into neuro-evolutionary algorithms to refine the optimization of neural networks. This method combines evolutionary strategies with constraint-driven feedback mechanisms, ensuring the evolution process adheres to a set of dynamically imposed mathematical and computational constraints. The evolution of neural network architectures and weights follows a rigorously constrained adaptation process, driven by feedback mechanisms that enforce both explicit and implicit constraints. Neuro-evolution optimizes the

Example ICONE Genome Visualization

<figure><img src="image_221.png" /><figcaption>Figure 16.52: Visualization of an Interactively Constrained Neuro-Evolution (ICONE) genome: Neurons are colored by type (input, hidden, output). Frozen neurons (under user constraint) are shaded gray. Arrows represent weighted connections (blue for positive, red for negative), with width proportional to weight magnitude. Disabled connections appear faded, and constrained connections are marked with an “X”</figcaption></figure>

weights and structures of artificial neural networks using evolutionary algorithms. Consider a neural network parameterized by a weight vector $\mathbf{w} = (w_1, w_2, \dots, w_n)$, where $n$ is the number of trainable parameters. The loss function of the network, denoted as $\mathcal{L}(\mathbf{w})$, defines the objective landscape in which evolutionary search occurs. The ICONE method introduces interactive constraints that enforce adherence to specific functional, structural, or performance-based criteria, which we define mathematically as

$$C_i(\mathbf{w}) \le 0, \quad i = 1, 2, \dots, m, \qquad (16.430)$$

where $C_i(\mathbf{w})$ represents the $i$th constraint function, ensuring the evolved network satisfies predefined conditions. The evolution process follows a mutation-selection cycle under constrained optimization. Given an initial population $\mathcal{P}^{(t)} = \{\mathbf{w}_1^{(t)}, \mathbf{w}_2^{(t)}, \dots, \mathbf{w}_N^{(t)}\}$ at generation $t$, candidate solutions undergo mutation:

$$\mathbf{w}_j^{(t+1)} = \mathbf{w}_j^{(t)} + \eta_j, \quad \eta_j \sim \mathcal{N}(0, \sigma^2 I), \qquad (16.431)$$
where $\eta_j$ is a perturbation vector sampled from an isotropic Gaussian distribution with variance $\sigma^2$. The feasibility of each mutated candidate is determined by evaluating the constraints $C_i(\mathbf{w}_j^{(t+1)})$. If a candidate violates any constraints, a projection step enforces feasibility by solving the constrained optimization problem:

$$ \mathbf{w}_j^{(t+1)} = \arg \min_{\mathbf{w}} \|\mathbf{w} - \mathbf{w}_j^{(t+1)}\|^2 \quad \text{s.t.} \quad C_i(\mathbf{w}) \le 0, \quad i = 1, 2, \dots, m. \qquad (16.432) $$

This ensures that every evolved individual remains within the feasible search space. Selection in ICONE is governed by fitness evaluation and constraint satisfaction. The fitness function $F(\mathbf{w})$ incorporates both performance metrics (e.g., classification accuracy, regression error) and constraint penalties. A penalty-based fitness function is formulated as:

$$ F(\mathbf{w}) = \mathcal{L}(\mathbf{w}) + \lambda \sum_{i=1}^{m} \max(0, C_i(\mathbf{w}))^2, \qquad (16.433) $$

where $\lambda$ is a penalty coefficient enforcing constraint satisfaction. Higher fitness candidates are selected for reproduction, forming the next generation $\mathcal{P}^{(t+1)}$. A key characteristic of ICONE is interactive constraint adaptation, where constraints evolve dynamically based on intermediate feedback. If the optimization process trends toward infeasible solutions, adaptive constraints are imposed by modifying the constraint functions:

$$ C_i^{(t+1)}(\mathbf{w}) = C_i^{(t)}(\mathbf{w}) + \gamma \cdot \nabla_{\mathbf{w}} C_i(\mathbf{w}), \qquad (16.434) $$

where $\gamma$ is an adaptation rate controlling the magnitude of constraint adjustment. This adaptation mechanism ensures the evolutionary process remains both stable and effective, guiding solutions toward desirable regions in the search space. The convergence of ICONE relies on satisfying the Karush-Kuhn-Tucker (KKT) conditions, which characterize optimality in constrained optimization. At convergence, optimal solutions satisfy the Lagrangian formulation:

$$ \mathcal{L}(\mathbf{w}, \boldsymbol{\mu}) = F(\mathbf{w}) + \sum_{i=1}^{m} \mu_i C_i(\mathbf{w}), \qquad (16.435) $$

where $\mu_i \ge 0$ are the Lagrange multipliers. The necessary conditions for optimality are:

$$ \nabla_{\mathbf{w}} \mathcal{L}(\mathbf{w}, \boldsymbol{\mu}) = 0, \quad \mu_i C_i(\mathbf{w}) = 0, \quad C_i(\mathbf{w}) \le 0, \quad \mu_i \ge 0. \qquad (16.436) $$

These conditions ensure the final neural network configuration is both performance-optimized and constraint-compliant.

Ultimately, the ICONE method provides a mathematically rigorous framework for constrained neuro-evolution, ensuring neural networks evolve under explicit, dynamically adaptable, and interactive constraints while maintaining optimal performance in their designated tasks. Through a constraint-driven evolutionary search, ICONE guarantees that evolved models satisfy both functional and theoretical constraints, leading to robust and interpretable artificial intelligence systems.

### 16.5.9.1 Python Code to generate Figure 16.52

Interactively Constrained Neuro-Evolution (ICONE) is a neuroevolution framework where:

1. The user (human-in-the-loop) constrains evolution by enforcing rules (e.g., disallow certain connections, force modularity, or bias toward certain structures).

2. Architectures evolve under both evolutionary pressure and interactive constraints.

3. Visualization should therefore highlight:
(a) Neurons and their types.

(b) Connections (weights, enabled/disabled).

(c) Constraints applied (e.g., disallowed connections, frozen neurons).

So, in addition to the usual NEAT/SUNA/DXNN-style graph, we need to mark constrained elements (e.g., red crosses on forbidden edges, shaded neurons if frozen). The Python code below produces the Figure 16.52 illustrating Interactively Constrained Neuro-Evolution (ICONE).

```python
"""
Plotting helper for ICONE (Interactively Constrained Neuro-Evolution).

Dependencies:
    pi install networkx matplotlib

"""

import matplotlib.pyplot as plt
import networkx as nx

# Example ICONE genome
example_genome = {
    "neurons": [
        {"id": "n1", "type": "input", "activation": "linear"},
        {"id": "n2", "type": "input", "activation": "linear"},
        {"id": "n3", "type": "hidden", "activation": "tanh"},
        {"id": "n4", "type": "hidden", "activation": "sigmoid"},
        {"id": "n5", "type": "output", "activation": "tanh"},
    ],
    "connections": [
        {"src": "n1", "dst": "n3", "w": 0.7, "enabled": True, "constrained": False},
        {"src": "n2", "dst": "n3", "w": -0.5, "enabled": True, "constrained": True},  # forbidden
        {"src": "n3", "dst": "n4", "w": 1.0, "enabled": True, "constrained": False},
        {"src": "n4", "dst": "n5", "w": 0.8, "enabled": True, "constrained": False},
        {"src": "n1", "dst": "n5", "w": 0.3, "enabled": False, "constrained": False},  # disabled
    ],
    "frozen_neurons": ["n4"]  # neuron under user constraint
}

# Visualization function

def plot_icone(genome, title="ICONE Genome", figsize=(9, 6), savefile=None):
    G = nx.DiGraph()

    # Add neurons
    for n in genome["neurons"]:
        G.add_node(n["id"], **n)
```
```python
# Add connections
for c in genome["connections"]:
    G.add_edge(c["src"], c["dst"], **c)

# Layout: spring
pos = nx.spring_layout(G, seed=42)

# Node properties
type_to_color = {"input": "#1f78b4", "hidden": "#fb9a99", "output":
                 "#33a02c"}
node_colors = []
node_shapes = []

for nid, attr in G.nodes(data=True):
    c = type_to_color.get(attr["type"], "gray")
    if nid in genome.get("frozen_neurons", []):
        c = "#aaaaaa"  # frozen = gray
    node_colors.append(c)

# Draw nodes
nx.draw_networkx_nodes(G, pos, node_color=node_colors, node_size
                      =1200,
                      edgecolors="k", linewidths=1.2)

# Draw labels
labels = {nid: f"{nid}\n{attr['activation']}" for nid, attr in G.
          nodes(data=True)}
nx.draw_networkxs(G, pos, labels=labels, font_size=9)

# Draw edges
ax = plt.gca()

for u, v, data in GEdges(data=True):
    w = data["w"]
    enabled = data["enabled"]
    constrained = data.get("constrained", False)

    color = "blue" if w >= 0 else "red"
    alpha = 1.0 if enabled else 0.2
    width = 1.0 + 2.0 * abs(w)

    ax annotate("",
              xy(pos[v], xycoordinates="data",
                  xytext(pos[u], textcoordinates="data",
                          arrowprops=dict(arrowstyle="->",
                                          linewidth=width,
                                          color=color,
                                          alpha=alpha,
                                          shrinkA=12, shrinkB=12)))

# If constrained: draw an "X" on the edge midpoint
if constrained:
    mx = (pos[u][0] + pos[v][0]) / 2
    my = (pos[u][1] + pos[v][1]) / 2
    ax.text(mx, my, "X", color="black", fontsize=12, ha="center", va="center")
```
<figure><img src="image_222.png" /></figure>

16.5.10 Deus Ex Neural Network (DXNN)

The Deus Ex Neural Network (DXNN) method is a sophisticated approach in artificial intelligence that employs a memetic algorithm-based Topology and Weight Evolving Artificial Neural Network (TWEANN) system. This method is designed to evolve both the structure and parameters of neural networks, enabling the development of models that can adapt to complex tasks without predefined architectures. In traditional neural network training, the architecture is often fixed, and only the

<figure><img src="image_223.png" /><figcaption>Figure 16.53: Visualization of a Deus Ex Neural Network (DXNN) genome. Neurons are grouped into cores (dashed rectangles), colored by type (input, hidden, output). Arrows represent weighted connections (blue for positive, red for negative), with line thickness proportional to weight magnitude. Disabled connections appear faded. Each neuron label indicates its identifier and activation function.</figcaption></figure>
weights are optimized using algorithms like backpropagation. In contrast, DXNN simultaneously evolves the topology and weights of the network. This is achieved through a combination of evolutionary strategies and local search methods, which iteratively refine the network's performance. The memetic algorithm aspect of DXNN incorporates local optimization techniques within the evolutionary process, enhancing convergence rates and solution quality. Mathematically, the DXNN method can be described as follows:

1. **Initialization**: A population of neural networks is initialized with random topologies and weights. Each network's topology can be represented as a graph $G = (V, E)$, where $V$ denotes neurons and $E$ denotes synaptic connections.

2. **Fitness Evaluation**: Each network $i$ in the population is evaluated based on a fitness function $F_i$, which measures its performance on a given task. The fitness function could be defined as:

$$F_i = \frac{1}{N} \sum_{j=1}^{N} L(y_j, \hat{y}_j) \qquad (16.437)$$

where $N$ is the number of samples, $y_j$ is the true output, $\hat{y}_j$ is the network's output, and $L$ is a loss function, such as mean squared error:

$$L(y, \hat{y}) = (y - \hat{y})^2 \qquad (16.438)$$

3. **Selection**: Networks are selected for reproduction based on their fitness scores. A common selection method is tournament selection, where a subset of networks is chosen, and the one with the highest fitness is selected for reproduction.

4. **Crossover (Recombination)**: Pairs of selected networks undergo crossover to produce offspring. This involves combining the topologies and weights of parent networks. For example, given parent networks with weight matrices $W^A$ and $W^B$, an offspring's weight matrix $W^O$ could be:

$$W^O = \alpha W^A + (1 - \alpha) W^B \qquad (16.439)$$

where $\alpha$ is a crossover coefficient.

5. **Mutation**: Offspring networks undergo mutations to introduce variability. Mutations can affect both the topology and weights. For weight mutation:

$$W' = W + \Delta W \qquad (16.440)$$

where $\Delta W$ is a perturbation matrix, often sampled from a normal distribution:

$$\Delta W \sim \mathcal{N}(0, \sigma^2) \qquad (16.441)$$

For topology mutation, connections can be added or removed based on a probability $p_{add}$ or $p_{remove}$.

6. **Local Optimization (Memetic Component)**: After mutation, local search methods, such as gradient-based optimization, are applied to fine-tune the weights of the offspring networks. This involves minimizing the loss function $L$ with respect to the weights $W$:
$$W \leftarrow W - \eta \nabla_{\theta} L \qquad (16.442)$$

where $\eta$ is the learning rate, and $\nabla_{\theta} L$ is the gradient of the loss function with respect to the weights.

7. **Replacement**: The new generation of networks replaces the old population, and the process repeats from the fitness evaluation step until a termination criterion is met, such as a predefined number of generations or a satisfactory fitness level.

The DXNN method has been applied in various domains, including financial markets. For instance, in automated currency trading, DXNN has been used to evolve neural networks that process Forex chart images as inputs, enabling the detection of patterns and trends for trading decisions. This approach contrasts with traditional methods that rely on fixed technical indicators, offering a more dynamic and adaptive trading strategy.

In summary, the Deus Ex Neural Network method provides a comprehensive framework for evolving both the architecture and parameters of neural networks. By integrating evolutionary algorithms with local optimization techniques, DXNN facilitates the development of adaptable and efficient models capable of tackling complex tasks across various domains.

### 16.5.10.1 Python Code to generate Figure 16.53

The Deus Ex Neural Network (DXNN) framework is a modular, hierarchical, evolving neural architecture. It differs from NEAT/SUNA in that: It is hierarchical:

1. Neurons form Cores (like sub-networks).

2. Cores are grouped into a DXNN organism.

Both neurons and cores can mutate, evolve, and rewire. Visualization should therefore show two levels:

1. The Core graph (inter-core connectivity).

2. Inside each core, the Neuron graph (intra-core connectivity).

3. Cores are grouped into a DXNN organism.

The Python code below produces the Figure 16.53 illustrating Deus Ex Neural Network (DXNN).

```python
"""
Plotting helper for Deus Ex Neural Network (DXNN).

Dependencies:
    pi install networkx matplotlib

This script visualizes a DXNN genome in a hierarchical manner:
- Cores are shown as big clusters (rectangles/ellipses).
- Neurons are shown inside their core cluster.
- Core-to-core connections and intra-core neuron connections are drawn.

Example genome is provided at the end.
"""
import matplotlib.pyplot as plt
```
```python
for n in genome["neurons"]:
    cores[n["core"]).append(n["id"])

pos = {}
core positions = {}
n cores = len(cores)
for i, (cid, nids) in enumerate(cores.items()):
    # Position cores along x-axis
cx = i * 4.0
cy = 0
core positions[cid] = (cx, cy)
# Arrange neurons vertically inside each core
for j, nid in enumerate(nids):
    pos[nid] = (cx, cy - j * 1.5)

# Draw neurons
type_to_color = {
    "input": "#1f78b4",
    "hidden": "#fb9a99",
    "output": "#33a02c",
}
node_colors = []
for n in G.nodes(data=True):
    node_colors.append(type_to_color.get(n[1]["type"], "gray"))
nx.draw_networkx_nodes(G, pos, node_color=node_colors, node_size=900,
                       edgecolors="k", linewidths=1.0, ax=ax)

# Draw labels (id + activation)
labels = {nid: f"{nid}\n{attr('activation')}" for nid, attr in G.
          nodes(data=True)}
nx.draw_networkx_labels(G, pos, labels=labels, font_size=8, ax=ax)

# Draw edges
for u, v, data in GEdges(data=True):
    w = data["weight"]
    enabled = data["enabled"]
    color = "blue" if w >= 0 else "red"
    alpha = 1.0 if enabled else 0.2
    width = 1.0 + 2.0 * abs(w)
    ax annotate("",
              xy = pos[v], xy coordinates='data',
              xytext = pos[u], text coordinates='data',
              arrowprops = dict(arrowstyle="->",
                               linewidth = width,
                               color = color,
                               alpha = alpha,
                               shrink A = 12,
                               shrink B = 12)))

# Draw core boundaries
for cid, (cx, cy) in core positions.items():
    nodes = cores[cid]
    if not nodes:
```
```python
continue
xs = [pos[n][0] for n in nodes]
ys = [pos[n][1] for n in nodes]
min_x, max_x = min(xs) - 1.0, max(xs) + 1.0
min_y, max_y = min(ys) - 1.0, max(ys) + 1.0
rect = mpatches parallelogram((min_x, min_y),
                             (max_x - min_x, max_y - min_y),
                             fill=False, linewidth=2,
                             edgecolor="black", alpha=0.5)
ax.add_patch(rect)
ax.text(cx, max_y + 0.5, cid, fontweight='bold', ha='center', va='bottom', fontweight='bold')

plt.title('title')
plt.axis('off')
plt.tight_layout()
plt.show()

# Example run
if __name__ == "__main__":
    plot_txnn(example_genome, title="Example DXNN Genome Visualization")
```

### 16.5.11 Spectrum-diverse Unified Neuroevolution Architecture (SUNA)

The Spectrum-diverse Unified Neuroevolution Architecture (SUNA) is an advanced framework in artificial intelligence that synergistically combines a unified neural representation with a novel diversity-preserving mechanism known as spectrum diversity. This integration facilitates the evolution of neural networks with both adaptable topologies and weights, enabling the system to proficiently address a wide array of complex tasks.

Central to SUNA is the Unified Neuron Model, which encapsulates various neuron types and activation functions within a single, cohesive representation. Each neuron *i* is characterized by a set of parameters $\theta_i$, governing its specific function. The output $y_i$ of neuron *i* can be mathematically expressed as:

$$y_i = f_i \left( \sum_{j \in \text{pre}(i)} w_{ij} x_j + b_i \right) \qquad (16.443)$$

Here, $f_i$ denotes the activation function, $w_{ij}$ represents the synaptic weight between neurons *j* and *i*, $x_j$ signifies the input from neuron *j*, and $b_i$ is the bias term associated with neuron *i*. This formulation allows for the seamless integration of diverse neuronal behaviors within a unified framework. The evolutionary process in SUNA involves the optimization of both the neural network's architecture and its synaptic weights. A population of candidate solutions, each encoded as a chromosome *C*, undergoes iterative selection, crossover, and mutation operations. The fitness $F(C)$ of each chromosome is evaluated based on its performance on the target task. The mutation operators are designed to modify the network's topology and weights, introducing variations that enhance the search for optimal solutions. To maintain a rich diversity of solutions, SUNA employs the spectrum diversity mechanism. This approach constructs a spectrum $S(C)$ for each chromosome, capturing
<figure><img src="image_224.png" /><figcaption>Figure 16.54: Visualization of the Spectrum-diverse Unified Neuroevolution Architecture (SUNA) genome: Node colors represent neuron types, node sizes correspond to neuron speeds, and edge widths/colors denote weight magnitude and sign. Disabled connections are shown with faded edges, while control neurons are highlighted with dashed borders</figcaption></figure>

its unique characteristics. The spectrum is defined as:

$$S(C) = (s_1(C), s_2(C), \dots, s_n(C)) \qquad (16.444)$$

where $s_k(C)$ represents the k-th feature of the chromosome C. The distance D between two spectra $S(C_1)$ and $S(C_2)$ is computed using a suitable metric, such as the Euclidean distance:

$$D(S(C_1), S(C_2)) = \sqrt{\sum_{k=1}^{n} (s_k(C_1) - s_k(C_2))^2} \qquad (16.445)$$

This distance metric informs the niching mechanism, ensuring that the evolutionary process explores a diverse set of solutions by promoting chromosomes with unique spectra. The fitness evaluation in SUNA is augmented by a novelty score $N(C)$, which quantifies the distinctiveness of a chromosome relative to the current population. The novelty score is calculated as:

$$N(C) = \frac{1}{k} \sum_{i=1}^{k} D(S(C), S(C_i)) \qquad (16.446)$$

where $C_i$ denotes the i-th nearest neighbor to C in the spectrum space, and k is a predefined constant. This scoring system encourages the exploration of novel solutions, thereby enhancing the algorithm's ability to escape local optima. The overall selection probability $P(C)$ of a chromosome is influenced by both its fitness and novelty scores, and can be expressed as:

$$P(C) = \frac{\alpha F(C) + \beta N(C)}{\sum_{C'} (\alpha F(C') + \beta N(C'))} \qquad (16.447)$$

Here, $\alpha$ and $\beta$ are weighting factors that balance the contributions of fitness and novelty, respectively. This probabilistic selection mechanism ensures a harmonious trade-off between exploiting
high-performing solutions and exploring innovative ones.

Through the integration of the Unified Neuron Model and spectrum diversity, SUNA adeptly navigates the complex search space inherent in neuroevolution. This comprehensive approach enables the discovery of neural network configurations that are both diverse and well-suited to a multitude of tasks, thereby advancing the field of artificial intelligence.

### 16.5.11.1 Python Code to generate Figure 16.54

The Python script below produces a clear, publication-quality visualization of a SUNA (Spectrum-diverse Unified Neuroevolution Architecture) genome which is given Figure 16.54. I make one small assumption: the genome is represented as a Python dict with a list of neurons and a list of connections (an example genome is included). The plot highlights different neuron types, node "speed" as size, control neurons, connection sign/strength, and disabled connections.

```python
#"
Plotting helper for SUNA-like genomes (Spectrum-diverse Unified Neuroevolution Architecture).
```

```
Dependencies:
```

```
  πp install networkx matplotlib
```

```
Usage:
```

```
  - Modify the 'example_genome' below to match your genome representation.
```

```
  - Run: python plot_suna.py
```

```python
#"
```

```python
import math
```

```python
import networkx as nx
```

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

```python
# ------------------------
```

```python
# Example SUNA genome
```

```python
# ------------------------
```

```python
example_genome = {
```

```python
    "neurons": [
        # id, type, activation, speed, controlInputs (ids)
        {"id": 0, "type": "input", "activation": "linear", "speed":
            1.0, "control": []},
        {"id": 1, "type": "input", "activation": "linear", "speed":
            1.0, "control": []},
        {"id": 2, "type": "hidden", "activation": "tanh", "speed":
            0.7, "control": []},
        {"id": 3, "type": "hidden", "activation": "sigmoid", "speed":
            0.4, "control": []},
        {"id": 4, "type": "control", "activation": "step", "speed":
            1.0, "control": []},
        {"id": 5, "type": "memory", "activation": "identity", "speed":
            0.2, "control": []},
        {"id": 6, "type": "output", "activation": "tanh", "speed":
            1.0, "control": [4]},
    ],
```

```python
    "connections": [
        # src, dst, weight, enabled
```
```python
        {"src": 0, "dst": 2, "w": 0.8, "enabled": True},
        {"src": 1, "dst": 2, "w": -0.4, "enabled": True},
        {"src": 2, "dst": 3, "w": 0.3, "enabled": True},
        {"src": 3, "dst": 6, "w": 1.2, "enabled": True},
        {"src": 5, "dst": 6, "w": 0.6, "enabled": True},
        {"src": 4, "dst": 2, "w": 1.0, "enabled": True}, # control ->
        hidden
        {"src": 1, "dst": 5, "w": 0.5, "enabled": False}, # disabled
        conn (plotted faded)
    ]
}
```
```python
for i, nid in enumerate inputs):
    pos[nid] = (left_x, in_pos[nid][1])

for i, nid in enumerate outputs):
    pos[nid] = (right_x, out_pos[nid][1])

# spread mid nodes around middle_x with slight offset so they don't
overlap with inputs/outputs
mid_x = middle_x

for i, nid in enumerate others):
    # put memory and control slightly above/below for clarity
    y = oth_pos[nid][1]
    pos[nid] = (mid_x + (0.12 * (i % 2)), y)

# Node visual properties: color by type, size by speed
type_to_color = {
    "input": "#1f78b4",
    "output": "#33a02c",
    "hidden": "#fb9a99",
    "control": "#6a3d9a",
    "memory": "#ff7f00",
    "other": "#b2df8a",
}
node_colors = []
node_sizes = []
node_edge styles = []
labels = {}

for n in G.nodes(data=True):
    nid, attr = n
    t = attr.get("type", "other")
    c = type_to_color.get(t, type_to_color["other"])
    node_colors.append(c)
    # size: base scaled by speed attribute (mapped 0.2->0.6)
    speed = float(attr.get("speed", 1.0))
    size = 300 + 1200 * (max(0.01, min(1.0, speed))) # mmapped to a
    reasonable range
    node_sizes.append(size)

# dashed border for control neurons to emphasize them
if t == "control":
    node_edge styles.append("dashed")
else:
    node_edge styles.append("solid")

# concise node label showing id and activation
labels[nid] = f"{nid}\n{attr.get('activation', '')}[:6]}"

# Draw base nodes (we'll overlay dashed borders if needed)
plt.figure(figsize=figsize)
ax = plt.gca()
nx.draw_networkx_nodes(G, pos,
    node_color=node_colors,
    node_size=node_sizes,
    linewidths=1.2,
    edgecolors="k",
    ax=ax)
```
```python
# overlay dashed borders for control nodes (by drawing them again
with no fill)
# Draw base nodes (filled)
plt.figure(figsize=figsize)
ax = plt.gca()
nx.draw_networkx_nodes(G, pos,
                       node_color=node_colors,
                       node_size=node_sizes,
                       linewidths=1.2,
                       edgecolors="k",
                       ax=ax)
```

```python
# Overlay dashed borders for control neurons
import matplotlib.patches as mpatches
for nid, style in zip(G.nodes(), node_edge_styles):
    if style == "dashed":
        x, y = pos[nid]
        r = math sqrt(node_sizes [list(G.nodes()).index(nid)]) /
        200.0 # scale radius from node size
        circ = mpatches-circle((x, y), radius=r,
                                 fill=False, linewidth="--", linewidth
                                 =2, edgecolor="k")
        ax.add_patch(circ)
```

```python
# Draw labels
nx.draw_networkx_labels(G, pos, labels=labels, font_size=9, ax=ax)
```

```python
# Draw edges, with color by sign and width by magnitude. Disabled
edges faded.
for u, v, data in GEdges(data=True):
    w = data.get("weight", 0.0)
    enabled = data.get("enabled", True)
    sign = "positive" if w >= 0 else "negative"
    mag = abs(w)
    # Normalize width within a good visual range
    width = 0.5 + 4.5 * (1 - math.exp(-abs(mag))) # sigmoidish
    scaling
    alpha = 1.0 if enabled else 0.18
    color = (0.8, 0.0, 0.0) if w < 0 else (0.0, 0.2, 0.6)
    # draw with arrow
    ax annotate("",
              xy=pos[v], xycoordinates='data',
              xytext=pos[u], textcoordinates='data',
              arrowprops=dict(arrowstyle="->",
                              linewidth=width,
                              color=color,
                              alpha=alpha,
                              shrinkA=10,
                              shrinkB=10),
              )
```

```python
# Legend creation
from matplotlib.lines import Line2D
legend_handles = []
```
```python
184 for t, c in type_to_color.items():
185     legend_handles.append(Line2D([0], [0], marker='o', color='w',
186                             label=t,
187                             markerfacecolor=c, markersize=10,
188                             markeredgecolor='k'))
187 # custom handle for disabled
188 legend_handles.append(Line2D([0], [0], linestyle='-', color
189                         =(0.5,0.5,0.5), label='disabled conn', linewidth=3, alpha=0.4))
190 plt.legend(handles=legend_handles, bbox_to塞 anchor=(1.02, 1.0), loc='upper left')
191 plt.title(title)
192 plt.axis('off')
193 plt.tight_layout()
194 plt.show()
195
196
197 # ----------------------------
198 # Example run
199 # ----------------------------
200 if __name__ == "__main__":
201     plot_suna(example_genome, title="Example SUNA Genome Visualization")
```
# 17 Training Neural Networks

Training neural networks involves optimizing a parameterized function $f_{\theta}: \mathcal{X} \to \mathcal{Y}$ to minimize an expected loss $\mathcal{L}(\theta) = \mathbb{E}_{(x,y)\sim\mathcal{D}}[\ell(f_{\theta}(x), y)]$, where $\mathcal{D}$ is the data distribution, $\ell$ is a per-sample loss function (e.g., cross-entropy or mean squared error), and $\theta \in \mathbb{R}^p$ represents the network's parameters. The optimization is typically performed via gradient-based methods, where the parameters are iteratively updated as:

$$ \theta_{t+1} = \theta_t - \eta_t \nabla_{\theta} \mathcal{L}(\theta_t), $$

with $\eta_t$ denoting the learning rate at iteration $t$. The gradient $\nabla_{\theta} \mathcal{L}$ is estimated using backpropagation, which decomposes the gradient computation via the chain rule. For a network with $L$ layers, the gradient of the loss with respect to the weights $\mathbf{W}^{(l)}$ of layer $l$ is given by:

$$ \nabla_{\mathbf{W}^{(l)}} \mathcal{L} = \mathbf{h}^{(l-1)} (\nabla_{\mathbf{z}^{(l)}} \mathcal{L})^T, $$

where $\mathbf{h}^{(l-1)}$ is the activation of the previous layer, $\mathbf{z}^{(l)} = \mathbf{W}^{(l)}\mathbf{h}^{(l-1)} + \mathbf{b}^{(l)}$ is the pre-activation, and $\nabla_{\mathbf{z}^{(l)}} \mathcal{L}$ is the error signal propagated backward from subsequent layers. The error signal for layer $l$ is computed recursively as:

$$ \nabla_{\mathbf{z}^{(l)}} \mathcal{L} = (\mathbf{W}^{(l+1)})^T \nabla_{\mathbf{z}^{(l+1)}} \mathcal{L} \odot \sigma'(\mathbf{z}^{(l)}), $$

where $\odot$ denotes element-wise multiplication and $\sigma'$ is the derivative of the activation function.

The training dynamics of neural networks are governed by the geometry of the loss landscape $\mathcal{L}(\theta)$, which is non-convex and high-dimensional. Stochastic gradient descent (SGD) introduces noise due to mini-batch sampling, leading to the parameter update:

$$ \theta_{t+1} = \theta_t - \eta_t (\nabla_{\theta} \mathcal{L}(\theta_t) + \xi_t), $$

where $\xi_t$ is a zero-mean noise term with covariance $\mathbf{C}(\theta_t)$. This noise can help escape saddle points and find flatter minima, which are empirically associated with better generalization. The continuous-time limit of SGD is described by the stochastic differential equation:

$$ d\theta(t) = -\nabla_{\theta} \mathcal{L}(\theta(t)) dt + \sqrt{\eta \mathbf{C}(\theta(t))} d\mathbf{W}(t), $$

where $\mathbf{W}(t)$ is a Wiener process. The diffusion term $\sqrt{\eta \mathbf{C}}$ influences the exploration of the parameter space, with larger learning rates $\eta$ amplifying the noise.

Regularization techniques are employed to prevent overfitting and improve generalization. Weight decay, equivalent to $L_2$-regularization, modifies the loss function as:

$$ \tilde{\mathcal{L}}(\theta) = \mathcal{L}(\theta) + \frac{\lambda}{2} \|\theta\|_2^2, $$

where $\lambda$ controls the strength of regularization. Dropout, another common technique, randomly masks activations during training, effectively approximating an ensemble of sub-networks. At test time, dropout is disabled, and the weights are scaled to account for the missing units. The implicit
bias of gradient descent refers to the tendency of optimization algorithms to converge to specific solutions even without explicit regularization. For example, in linear models, gradient descent converges to the minimum-norm solution when initialized at zero, while in nonlinear neural networks, the implicit bias often favors solutions with low complexity, as measured by the norm or margin.

The convergence of gradient-based optimization is analyzed using Lyapunov stability theory. Under smoothness and convexity assumptions, the loss decreases monotonically as:

$$ \mathcal{L}(\theta_{t+1}) \leq \mathcal{L}(\theta_t) - \frac{\eta}{2} \|\nabla_{\theta}\mathcal{L}(\theta_t)\|^2 + O(\eta^2). $$

In non-convex settings, gradient descent converges to critical points where $\nabla_{\theta}\mathcal{L}(\theta) = 0$. The nature of these critical points (minima, maxima, or saddle points) is determined by the Hessian $\nabla^2\mathcal{L}(\theta)$. Second-order methods, such as Newton's method, leverage the Hessian to accelerate convergence:

$$ \theta_{t+1} = \theta_t - \eta (\nabla^2\mathcal{L}(\theta_t))^{-1} \nabla_{\theta}\mathcal{L}(\theta_t), $$

but are computationally expensive for large networks. Approximate second-order methods, like Adam and RMSprop, adapt the learning rate per parameter, combining momentum and curvature information to improve training efficiency. The interplay between optimization dynamics, architecture design, and regularization defines the success of neural network training in practice.

## 17.1 Literature Review of Training Neural Networks

Sorrenson (2025) [128] introduced a framework enabling exact maximum likelihood training of unrestricted neural networks. It presents new training methodologies based on probabilistic models and applies them to scientific applications. Liu and Shi (2015) [129] applied advanced neural network theory to meteorological predictions. It uses sensitivity analysis and new training techniques to mitigate sample size limitations. Das et. al. (2025) [130] integrated Finite Integral Transform (FIT) with gradient-enhanced physics-informed neural networks (g-PINN), optimizing training in engineering applications. Zhang et. al. (2025) [131] in his thesis explores neural tangent kernel (NTK) theory to model the gradient descent training process of deep networks and its implications for structural identification. Ali and Hussein (2025) [132] developed a hybrid approach combining fuzzy set theory and artificial neural networks, enhancing training robustness through heuristic optimization. Li (2025) [133] introduced a deep learning-based strategy to train neural networks for imperfect-information extensive-form games, emphasizing offline training techniques. Ghosh (2025) [1066] presented a mathematically rigorous continuum plasticity framework enriched with gradient terms to capture nonlocal effects in beam deformation. The model remains rate-independent, preserving classical yield and hardening behavior, while the inclusion of gradient-enhanced variables facilitates the representation of size-dependent behavior and mitigates localization issues commonly encountered in standard plasticity theories. By formulating the governing equations specifically for Euler-Bernoulli beam kinematics, the study delivers a tractable yet robust methodology for problems where gradient effects play a critical role in plastic deformation. Ghosh (2025) [992] developed a thermodynamically consistent and geometrically exact constitutive framework for rate-dependent plasticity within the scope of Euler-Bernoulli beam theory. The model rigorously integrates rate-dependent plastic flow into beam kinematics, thereby enabling accurate representation of size effects and time-dependent plastic deformation under large strains and elevated temperatures. Comprehensive derivations of the constitutive relations, alongside numerical illustrations, showcase the robustness of the proposed model for capturing complex creep-like responses in slender structures. Hu et. al. (2025) [135] explored the convergence properties of deep learning-based PDE solvers, analyzing training loss and function space properties. Chen et. al. (2025) [136] developed a Transformer-based neural network training framework for risk analysis, incorporating feature maps and game-theoretic interpretation. Sun et. al. (2025) [137] established
a new benchmarking suite for optimizing neural architecture search (NAS) techniques in training spiking neural networks. Zhang et. al. (2025) [138] proposed a novel iterative training approach for neural networks, enhancing convergence guarantees in theory and practice.

## 17.2 Backpropagation Algorithm

Consider a neural network with $L$ layers, where each layer $l$ (with $l = 1, 2, \dots, L$) consists of a weight matrix $\mathbf{W}^{(l)} \in \mathbb{R}^{n_l \times n_{l-1}}$, a bias vector $\mathbf{b}^{(l)} \in \mathbb{R}^{n_l}$, and an activation function $\sigma^{(l)}$ which is applied element-wise. The network takes as input a vector $x^{(i)} \in \mathbb{R}^{n_0}$ for the $i$-th training sample, where $n_0$ is the number of input features, and propagates it through the layers to produce an output $\hat{y}^{(i)} \in \mathbb{R}^{n_L}$, where $n_L$ is the number of output units. The network parameters (weights and biases) $\theta = \{\mathbf{W}^{(l)}, \mathbf{b}^{(l)}\}_{l=1}^L$ are to be optimized to minimize a loss function that captures the error between the predicted output $\hat{y}^{(i)}$ and the true target $y^{(i)}$ for all training examples. For each training sample, we define the loss function $\mathcal{L}(\hat{y}^{(i)}, y^{(i)})$ as the squared error:

$$ \mathcal{L}(\hat{y}^{(i)}, y^{(i)}) = \frac{1}{2} \|\hat{y}^{(i)} - y^{(i)}\|_2^2, \qquad (17.1) $$

where $\|\cdot\|_2$ represents the Euclidean norm. The total loss $J(\theta)$ for the entire dataset is the average of the individual losses:

$$ J(\theta) = \frac{1}{N} \sum_{i=1}^{N} \mathcal{L}(\hat{y}^{(i)}, y^{(i)}), \qquad (17.2) $$

where $N$ is the number of training samples. For squared error loss, we can write:

$$ J(\theta) = \frac{1}{2N} \sum_{i=1}^{N} \|\hat{y}^{(i)} - y^{(i)}\|_{2}^{2}. \qquad (17.3) $$

The forward pass through the network consists of computing the activations at each layer. For the $l$-th layer, the pre-activation $z^{(l)}$ is calculated as:

$$ z^{(l)} = \mathbf{W}^{(l)} a^{(l-1)} + \mathbf{b}^{(l)}, \qquad (17.4) $$

where $a^{(l-1)}$ is the activation from the previous layer and $\mathbf{W}^{(l)}$ is the weight matrix connecting the $(l-1)$-th layer to the $l$-th layer. The output of the layer, i.e., the activation $a^{(l)}$, is computed by applying the activation function $\sigma^{(l)}$ element-wise to $z^{(l)}$:

$$ a^{(l)} = \sigma^{(l)}(z^{(l)}). \qquad (17.5) $$

The final output of the network is given by the activation $a^{(L)}$ at the last layer, which is the predicted output $\hat{y}^{(i)}$:

$$ \hat{y}^{(i)} = a^{(L)}. \qquad (17.6) $$

The backpropagation algorithm computes the gradient of the loss function $J(\theta)$ with respect to each parameter (weights and biases). First, we compute the error at the output layer. Let $\delta^{(L)}$ represent the error at layer $L$. This is computed by taking the derivative of the loss function with respect to the activations at the output layer:

$$ \delta^{(L)} = \frac{\partial \mathcal{L}}{\partial a^{(L)}} \odot \sigma^{(L)'}(z^{(L)}), \qquad (17.7) $$

where $\odot$ denotes element-wise multiplication, and $\sigma^{(L)'}(z^{(L)})$ is the derivative of the activation function applied element-wise to $z^{(L)}$. For squared error loss, the derivative with respect to the activations is:

$$ \frac{\partial \mathcal{L}}{\partial a^{(L)}} = \hat{y}^{(i)} - y^{(i)} \qquad (17.8) $$
<figure><img src="image_225.png" /><figcaption>Figure 17.1: Visualization of the Backpropagation Algorithm: forward activations (black arrows) and backward error propagation (red dashed arrows)</figcaption></figure>

so the error term at the output layer is:

$$ \delta^{(L)} = (\hat{y}^{(i)} - y^{(i)}) \odot \sigma^{(L)'}(z^{(L)}) \quad (17.9) $$

To propagate the error backward through the network, we compute the errors at the hidden layers. For each hidden layer $l = L - 1, L - 2, \dots, 1$, the error $\delta^{(l)}$ is calculated by the chain rule:

$$ \delta^{(l)} = (\mathbf{W}^{(l+1)T} \delta^{(l+1)}) \odot \sigma^{(l)'}(z^{(l)}) \quad (17.10) $$

where $\mathbf{W}^{(l+1)T} \in \mathbb{R}^{n_{l+1} \times n_l}$ is the transpose of the weight matrix connecting layer $l$ to layer $l+1$. This equation uses the fact that the error at layer $l$ depends on the error at the next layer, modulated by the weights, and the derivative of the activation function at layer $l$. Once the errors $\delta^{(l)}$ are computed for all layers, we can compute the gradients of the loss function with respect to the parameters (weights and biases). The gradient of the loss with respect to the weights $\mathbf{W}^{(l)}$ is:

$$ \frac{\partial J(\theta)}{\partial \mathbf{W}^{(l)}} = \frac{1}{N} \sum_{i=1}^{N} \delta^{(l)}(a^{(l-1)})^{T} \quad (17.11) $$

The gradient of the loss with respect to the biases $\mathbf{b}^{(l)}$ is:

$$ \frac{\partial J(\theta)}{\partial \mathbf{b}^{(l)}} = \frac{1}{N} \sum_{i=1}^{N} \delta^{(l)} \quad (17.12) $$

After computing these gradients, we update the parameters using an optimization algorithm such as gradient descent. The weight update rule is:

$$ \mathbf{W}^{(l)} \leftarrow \mathbf{W}^{(l)} - \eta \frac{\partial J(\theta)}{\partial \mathbf{W}^{(l)}}, \quad (17.13) $$

and the bias update rule is:

$$ \mathbf{b}^{(l)} \leftarrow \mathbf{b}^{(l)} - \eta \frac{\partial J(\theta)}{\partial \mathbf{b}^{(l)}} \quad (17.14) $$
where η is the learning rate controlling the step size in the gradient descent update. This process of forward pass, backpropagation, and parameter update is repeated over multiple epochs, with each epoch consisting of a forward pass, a backward pass, and a parameter update, until the network converges to a local minimum of the loss function.

At each step of backpropagation, the chain rule is applied recursively to propagate the error backward through the network, adjusting each weight and bias to minimize the total loss. The derivative of the activation function $\sigma^{(l)}(z^{(l)})$ is critical, as it dictates how the error is modulated at each layer. Depending on the choice of activation function (e.g., ReLU, sigmoid, or tanh), the derivative will take different forms, and this choice has a direct impact on the learning dynamics and convergence rate of the network. Thus, backpropagation serves as the computational backbone of neural network training. By calculating the gradients of the loss function with respect to the network parameters through efficient error propagation, backpropagation allows the network to adjust its parameters iteratively, gradually minimizing the error and improving its performance across tasks. This process is mathematically rigorous, utilizing fundamental principles of calculus and optimization, ensuring that the neural network learns effectively from its training data.

### 17.2.0.1 Python Code to generate Figure 17.1

This code will draw a 3-layer neural network (input, hidden, output) and annotate the forward activations and backpropagation error signals.

```python
import matplotlib.pyplot as plt
import networkx as nx

# Define the neural network structure
layers = [3, 4, 2]  # 3 input neurons, 4 hidden neurons, 2 output neurons

# Create a directed graph
G = nx.DiGraph()

# Positioning of nodes
pos = {}
layer-gap = 3
node-gap = 1.5

for l, num_nodes in enumerate layers:
    for n in range(num_nodes):
        node_id = f"L{l}N{n}"
        pos[node_id] = (l * layer-gap, -n * node-gap)
        G.add_node(node_id, layer=l)

# Add forward edges (black)
for l in range(len layers) - 1:
    for n1 in range layers[l]:
        for n2 in range layers[l+1]:
            G.add_edge(f"L{l}N{n1}", f"L{l+1}N{n2}", color="black")

# Add backward edges (red, dashed)
for l in range(len layers) - 1, 0, -1:
    for n1 in range layers[l]:
        for n2 in range layers[l-1]:
            G.add_edge(f"L{l}N{n1}", f"L{l-1}N{n2}", color="red")
```
```python
33 # Draw forward edges
34 edge_colors = [G[u][v]>'color'] for u, v in GEdges]
35 styles = ['solid' if c == 'black' else 'dashed' for c in edge_colors]
36
37 plt.figure(figsize=(10,6))
38 for (u,v), style, color in zip(Gedges, styles, edge_colors):
39     nx.draw_networkx_edges(G, pos, edgelist=[(u,v)], edge_color=color,
40                           style=style, arrows=True, arrowsize=15)
41
42 # Draw nodes
43 nx.draw_networkx_nodes(G, pos, node_color="lightblue", node_size=1000)
44 nx.draw_networkx labels(G, pos, font_size=10, font_weight="bold")
45
46 plt.title("Backpropagation: Forward pass (black), Error backprop (red
dashed)")
47 plt axis("off")
48 plt show()
```

## 17.3 Gradient Descent Variants

The training of neural networks using gradient descent and its variants is a mathematically intensive process that aims to minimize a differentiable scalar loss function $\mathcal{L}(\theta)$, where $\theta$ represents the parameter vector of the neural network. The loss function is often expressed as

$$ \mathcal{L}(\theta) = \frac{1}{N} \sum_{i=1}^{N} \ell(\theta; x_i, y_i), \qquad (17.15) $$

where $(x_i, y_i)$ are the input-output pairs in the training dataset of size $N$, and $\ell(\theta; x_i, y_i)$ is the sample-specific loss. The minimization problem is solved iteratively, starting from an initial guess $\theta^{(0)}$ and updating according to the rule

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

where $\eta > 0$ is the learning rate, and $\nabla_{\theta} \mathcal{L}(\theta)$ is the gradient of the loss with respect to $\theta$. The gradient, computed via backpropagation, follows the chain rule and propagates through the network's layers to adjust weights and biases optimally. In a feedforward neural network with $L$ layers, the computations proceed as follows. The input to layer $l$ is

$$ z^{(l)} = W^{(l)} a^{(l-1)} + b^{(l)}, \qquad (17.17) $$

where $W^{(l)} \in \mathbb{R}^{n_l \times n_{l-1}}$ and $b^{(l)} \in \mathbb{R}^{n_l}$ are the weight matrix and bias vector for the layer, respectively, and $a^{(l-1)}$ is the activation vector from the previous layer. The output is then

$$ a^{(l)} = f^{(l)}(z^{(l)}), \qquad (17.18) $$

where $f^{(l)}$ is the activation function. Backpropagation begins with the computation of the error at the output layer,

$$ \delta^{(L)} = \frac{\partial \ell}{\partial a^{(L)}} \odot f'^{(L)}(z^{(L)}), \qquad (17.19) $$

where $f'^{(L)}(\cdot)$ is the derivative of the activation function. For hidden layers, the error propagates recursively as

$$ \delta^{(l)} = (W^{(l+1)})^{\top} \delta^{(l+1)} \odot f'^{(l)}(z^{(l)}). \qquad (17.20) $$
The gradients for weight and bias updates are then computed as

$$ \frac{\partial \mathcal{L}}{\partial \mathbf{W}^{(l)}} = \boldsymbol{\delta}^{(l)}(\boldsymbol{a}^{(l-1)})^{\top} \qquad (17.21) $$

and

$$ \frac{\partial \mathcal{L}}{\partial \mathbf{b}^{(l)}} = \boldsymbol{\delta}^{(l)}, \qquad (17.22) $$

respectively. The dynamics of gradient descent are deeply influenced by the curvature of the loss surface, encapsulated by the Hessian matrix

$$ \mathbf{H}(\boldsymbol{\theta}) = \nabla_{\boldsymbol{\theta}}^2 \mathcal{L}(\boldsymbol{\theta}). \qquad (17.23) $$

For a small step size $\eta$, the change in the loss function can be approximated as

$$ \Delta \mathcal{L} \approx -\eta \| \nabla_{\boldsymbol{\theta}} \mathcal{L}(\boldsymbol{\theta}) \|^2 + \frac{\eta^2}{2} (\nabla_{\boldsymbol{\theta}} \mathcal{L}(\boldsymbol{\theta}))^{\top} \mathbf{H}(\boldsymbol{\theta}) \nabla_{\boldsymbol{\theta}} \mathcal{L}(\boldsymbol{\theta}). \qquad (17.24) $$

This reveals that convergence is determined not only by the gradient magnitude but also by the curvature of the loss surface along the gradient direction. The eigenvalues $\lambda_1, \lambda_2, \dots, \lambda_d$ of $\mathbf{H}(\boldsymbol{\theta})$ dictate the local geometry, with large condition numbers $\kappa = \frac{\lambda_{\max}}{\lambda_{\min}}$ slowing convergence due to ill-conditioning. Stochastic gradient descent (SGD) modifies the standard gradient descent by computing updates based on a single data sample $(\boldsymbol{x}_i, y_i)$, leading to

$$ \boldsymbol{\theta}^{(k+1)} = \boldsymbol{\theta}^{(k)} - \eta \nabla_{\boldsymbol{\theta}} \ell(\boldsymbol{\theta}; \boldsymbol{x}_i, y_i). \qquad (17.25) $$

While SGD introduces variance into the updates, this stochasticity helps escape saddle points characterized by zero gradient but mixed curvature. To balance computational efficiency and stability, mini-batch SGD computes gradients over a randomly selected subset $\mathcal{B} \subset \{1, \dots, N\}$ of size $|\mathcal{B}|$, yielding

$$ \nabla_{\boldsymbol{\theta}} \mathcal{L}_{\mathcal{B}}(\boldsymbol{\theta}) = \frac{1}{|\mathcal{B}|} \sum_{i \in \mathcal{B}} \nabla_{\boldsymbol{\theta}} \ell(\boldsymbol{\theta}; \boldsymbol{x}_i, y_i). \qquad (17.26) $$

Momentum methods enhance convergence by incorporating a memory of past gradients. The velocity term

$$ \boldsymbol{v}^{(k+1)} = \gamma \boldsymbol{v}^{(k)} + \eta \nabla_{\boldsymbol{\theta}} \mathcal{L}(\boldsymbol{\theta}) \qquad (17.27) $$

accumulates gradient information, and the parameter update is

$$ \boldsymbol{\theta}^{(k+1)} = \boldsymbol{\theta}^{(k)} - \boldsymbol{v}^{(k+1)}. \qquad (17.28) $$

Analyzing momentum in the eigenspace of $\mathbf{H}(\boldsymbol{\theta})$, with $\mathbf{H} = \mathbf{Q}\mathbf{\Lambda}\mathbf{Q}^{\top}$, reveals that the effective step size in each eigendirection is

$$ \eta_{\text{eff},i} = \frac{\eta}{1 - \gamma \lambda_i}, \qquad (17.29) $$

showing that momentum accelerates convergence in low-curvature directions while damping oscillations in high-curvature directions. Adaptive gradient methods, such as AdaGrad, RMSProp, and Adam, refine learning rates for individual parameters. In AdaGrad, the adaptive learning rate is

$$ \eta_i^{(k+1)} = \frac{\eta}{\sqrt{G_{ii}^{(k+1)} + \epsilon}}, \qquad (17.30) $$

where

$$ \mathbf{G}_{ii}^{(k+1)} = \mathbf{G}_{ii}^{(k)} + (\nabla_{\boldsymbol{\theta}_i} \mathcal{L}(\boldsymbol{\theta}_i))^2. \qquad (17.31) $$

RMSProp modifies this with an exponentially weighted average

$$ \mathbf{G}_{ii}^{(k+1)} = \beta \mathbf{G}_{ii}^{(k)} + (1 - \beta) (\nabla_{\boldsymbol{\theta}_i} \mathcal{L}(\boldsymbol{\theta}_i))^2. \qquad (17.32) $$
Adam combines RMSProp with momentum, where the first and second moments are

$$ \mathbf{m}^{(k+1)} = \beta_1 \mathbf{m}^{(k)} + (1 - \beta_1) \nabla_{\boldsymbol{\theta}} \mathcal{L}(\boldsymbol{\theta}) \qquad (17.33) $$

and

$$ \mathbf{v}^{(k+1)} = \beta_2 \mathbf{v}^{(k)} + (1 - \beta_2) (\nabla_{\boldsymbol{\theta}} \mathcal{L}(\boldsymbol{\theta}))^2. \qquad (17.34) $$

Bias corrections yield

$$ \hat{\mathbf{m}}^{(k+1)} = \frac{\mathbf{m}^{(k+1)}}{1 - \beta_1^k}, \quad \hat{\mathbf{v}}^{(k+1)} = \frac{\mathbf{v}^{(k+1)}}{1 - \beta_2^k}. \qquad (17.35) $$

The final parameter update is

$$ \boldsymbol{\theta}^{(k+1)} = \boldsymbol{\theta}^{(k)} - \eta \frac{\hat{\mathbf{m}}^{(k+1)}}{\sqrt{\hat{\mathbf{v}}^{(k+1)} + \epsilon}}. \qquad (17.36) $$

In conclusion, gradient descent and its variants provide a rich framework for optimizing neural network parameters. While standard gradient descent offers a basic approach, advanced methods like momentum and adaptive gradients significantly enhance convergence by tailoring updates to the landscape of the loss surface and the dynamics of training.

## 17.3.1 SGD (Stochastic Gradient Descent) Optimizer

### 17.3.1.1 Literature Review of SGD (Stochastic Gradient Descent) Optimizer

Lauand and Meyn (2025) [186] established a theoretical framework for SGD using Markovian dynamics to improve convergence properties. It integrates quasi-periodic linear systems into SGD, enhancing its robustness in non-stationary environments. Maranjyan et al. (2025) [187] developed an asynchronous SGD algorithm that meets the theoretical lower bounds for time complexity. It introduces ring-based communication to optimize parallel execution without degrading convergence rates. Gao and Gündüz (2025) [188] proposed a stochastic gradient descent-based approach to optimize graph neural networks in wireless networks. It rigorously analyzes the stochastic optimization problem and proves its convergence guarantees. Yoon et al. (2025) [189] investigated federated SGD in multi-agent learning and derives theoretical guarantees on its communication efficiency while achieving equilibrium. Verma and Maiti (2025) [190] proposed a periodic learning rate (using sine and cosine functions) for SGD-based optimizers, theoretically proving its benefits in stability and computational efficiency. Borowski and Miasojedow (2025) [191] extended the Robbins-Monro theorem to analyze convergence guarantees of SGD, refining the theoretical understanding of projected stochastic approximation algorithms. Dong et al. (2025) [192] applied stochastic gradient descent to brain network modeling, providing a theoretical framework for optimizing neural control strategies. Jiang et al. (2025) [193] analyzed the bias-variance tradeoff in decentralized SGD, proving convergence rates and proposing an error-correction mechanism for biased gradients. Sonobe et al. (2025) [194] connected SGD with Bayesian inference, presenting a theoretical analysis of how stochastic optimization methods approximate posterior distributions. Ghosh (2020) [105] explored fundamental mathematical inequalities and their applications across different areas of analysis and algebra. It discusses various classical inequality forms, their proofs, and the conditions under which equality holds. The work aims to provide a clearer understanding of the relationships and bounds that govern numerical and functional expressions in mathematics. Zhang and Jia (2025) [195] examined the theoretical properties of policy gradients in reinforcement learning, proving convergence guarantees for stochastic optimal control problems.

### 17.3.1.2 Analysis of SGD (Stochastic Gradient Descent) Optimizer

The **Stochastic Gradient Descent (SGD) optimizer** is an iterative method designed to minimize an objective function $f(\mathbf{w})$ by updating a parameter vector $\mathbf{w}$ in the direction of the negative
gradient. The fundamental optimization problem can be expressed as

$$ \min_{\mathbf{w}} f(\mathbf{w}), \qquad (17.37) $$

where

$$ f(\mathbf{w}) = \frac{1}{N} \sum_{i=1}^{N} \ell(\mathbf{w}; \mathbf{x}_i, y_i) \qquad (17.38) $$

represents the empirical risk, constructed from a dataset $\{(\mathbf{x}_i, y_i)\}_{i=1}^N$. Here, $\ell(\mathbf{w}; \mathbf{x}_i, y_i)$ denotes the loss function, $\mathbf{w} \in \mathbb{R}^d$ is the parameter vector, $N$ is the dataset size, and $f(\mathbf{w})$ approximates the true population risk

$$ \mathbb{E}_{\mathbf{x},y}[\ell(\mathbf{w}; \mathbf{x}, y)]. \qquad (17.39) $$

Standard gradient descent involves the update rule

$$ \mathbf{w}^{(t+1)} = \mathbf{w}^{(t)} - \eta \nabla f(\mathbf{w}^{(t)}), \qquad (17.40) $$

where $\eta > 0$ is the learning rate and

$$ \nabla f(\mathbf{w}) = \frac{1}{N} \sum_{i=1}^{N} \nabla \ell(\mathbf{w}; \mathbf{x}_i, y_i) \qquad (17.41) $$

is the full gradient. However, for large-scale datasets, the computation of $\nabla f(\mathbf{w})$ becomes computationally prohibitive, motivating the adoption of stochastic approximations. The stochastic approximation relies on the idea of estimating the gradient $\nabla f(\mathbf{w})$ using a single data point or a small batch of data points. Denoting the random index sampled at iteration $t$ as $i_t$, the stochastic gradient can be written as

$$ \hat{\nabla} f(\mathbf{w}^{(t)}) = \nabla \ell(\mathbf{w}^{(t)}; \mathbf{x}_{i_t}, y_{i_t}). \qquad (17.42) $$

Consequently, the update rule becomes

$$ \mathbf{w}^{(t+1)} = \mathbf{w}^{(t)} - \eta \hat{\nabla} f(\mathbf{w}^{(t)}). \qquad (17.43) $$

For a mini-batch $\mathcal{B}_t$ of size $m$, the stochastic gradient generalizes to

$$ \hat{\nabla} f(\mathbf{w}^{(t)}) = \frac{1}{m} \sum_{i \in \mathcal{B}_t} \nabla \ell(\mathbf{w}^{(t)}; \mathbf{x}_i, y_i). \qquad (17.44) $$

An important property of $\hat{\nabla} f(\mathbf{w})$ is its unbiasedness:

$$ \mathbb{E}[\hat{\nabla} f(\mathbf{w})] = \nabla f(\mathbf{w}). \qquad (17.45) $$

However, the variance of $\hat{\nabla} f(\mathbf{w})$, defined as

$$ \mathrm{Var}[\hat{\nabla} f(\mathbf{w})] = \mathbb{E}[\|\hat{\nabla} f(\mathbf{w}) - \nabla f(\mathbf{w})\|^2], \qquad (17.46) $$

introduces stochastic noise into the updates, where $\mathrm{Var}[\hat{\nabla} f(\mathbf{w})] \approx \frac{\sigma^2}{m}$ and

$$ \sigma^2 = \mathbb{E}[\|\nabla \ell(\mathbf{w}; \mathbf{x}) - \nabla f(\mathbf{w})\|^2] \qquad (17.47) $$

is the variance of the gradients. To analyze the convergence properties of SGD, we assume $f(\mathbf{w})$ to be L-smooth, meaning

$$ \|\nabla f(\mathbf{w}_1) - \nabla f(\mathbf{w}_2)\| \le L \|\mathbf{w}_1 - \mathbf{w}_2\|, \qquad (17.48) $$

and $f(\mathbf{w})$ to be bounded below by $f^* = \inf_{\mathbf{w}} f(\mathbf{w})$. Using Taylor expansion, we can write

$$ f(\mathbf{w}^{(t+1)}) \le f(\mathbf{w}^{(t)}) - \eta \|\nabla f(\mathbf{w}^{(t)})\|^2 + \frac{\eta^2 L}{2} \|\hat{\nabla} f(\mathbf{w}^{(t)})\|^2. \qquad (17.49) $$
<figure><img src="image_226.png" /><figcaption>Figure 17.2: Trajectory of the Stochastic Gradient Descent (SGD) optimizer on $f(x, y) = x^2 + y^2$</figcaption></figure>

Taking expectations yields

$$ \mathbb{E}[f(\mathbf{w}^{(t+1)})] \le \mathbb{E}[f(\mathbf{w}^{(t)})] - \frac{\eta}{2}\mathbb{E}[\|\nabla f(\mathbf{w}^{(t)})\|^2] + \frac{\eta^2 L}{2}\sigma^2, \quad (17.50) $$

showing that the convergence rate depends on the interplay between the learning rate $\eta$, the smoothness constant $L$, and the gradient variance $\sigma^2$. For $\eta$ small enough, the dominant term in convergence is $-\frac{\eta}{2}\mathbb{E}[\|\nabla f(\mathbf{w}^{(t)})\|^2]$, leading to monotonic decrease in $f(\mathbf{w}^{(t)})$. In the strongly convex case, where $f(\mathbf{w})$ satisfies

$$ f(\mathbf{w}_1) \ge f(\mathbf{w}_2) + \nabla f(\mathbf{w}_2)^{\top}(\mathbf{w}_1 - \mathbf{w}_2) + \frac{\mu}{2}\|\mathbf{w}_1 - \mathbf{w}_2\|^2 \quad (17.51) $$

for $\mu > 0$, SGD achieves linear convergence. Specifically,

$$ \mathbb{E}[\|\mathbf{w}^{(t)} - \mathbf{w}^*\|^2] \le (1 - \eta\mu)^t \|\mathbf{w}^{(0)} - \mathbf{w}^*\|^2 + \frac{\eta\sigma^2}{2\mu}. \quad (17.52) $$

For non-convex functions, where $\nabla^2 f(\mathbf{w})$ can have both positive and negative eigenvalues, SGD may converge to a local minimizer or saddle point. Stochasticity plays a pivotal role in escaping strict saddle points $\mathbf{w}_s$ where $\nabla f(\mathbf{w}_s) = 0$ but $\lambda_{\min}(\nabla^2 f(\mathbf{w}_s)) < 0$.

### 17.3.1.3 Python Code to generate Figure 17.2

We'll again use the simple quadratic test function

$$ f(x, y) = x^2 + y^2 \quad (17.53) $$

to visualize the optimization trajectory.
import numpy as np
import matplotlib.pyplot as plt

# Objective function
def f(x, y):
    return x**2 + y**2

# Gradient of the function
def grad_f(x, y):
    return np.array([2*x, 2*y])

# SGD parameters
lr = 0.1  # learning rate
iterations = 20
x, y = 3.0, 3.0  # starting point

# Store the trajectory
trajectory = [(x, y)]

# SGD updates
for i in range(iterations):
    gradient = grad_f(x, y)
    # Adding a small stochastic noise
    noise = np.random.normal(0, 0.1, size=2)
    x, y = np.array([x, y]) - lr * gradient + noise
    trajectory.append((x, y))

trajectory = np.array(trajectory)

# Plotting
X = np.linspace(-4, 4, 100)
Y = np.linspace(-4, 4, 100)
X, Y = np.meshgrid(X, Y)
Z = f(X, Y)

plt.figure(figsize=(8,6))
plt.contour(X, Y, Z, levels=30, cmap='viridis')
plt plot(trajectory[:,0], trajectory[:,1], 'ro-', label='SGD Path')

# SGD Optimizer Trajectory on $f(x, y) = x^2 + y^2$

plt.title("SGD Optimizer Trajectory on $f(x, y) = x^2 + y^2$")

plt.tight_layout()

17.3.2 Nesterov Accelerated Gradient Descent (NAG)

17.3.2.1 Literature Review of Nesterov Accelerated Gradient Descent (NAG)

The field of Nesterov Accelerated Gradient Descent (NAG) has undergone significant theoretical
refinement and practical adaptation in recent years, with researchers delving into its convergence
properties, dynamical systems interpretations, stochastic extensions, and domain-specific optimiza-
tions. Adly and Attouch (2024) [440] provide an in-depth complexity analysis by precisely tuning
the viscosity parameter within an inertial gradient system, thereby extending NAG's classical for-
mulations into the Su-Boyd-Candès dynamical framework. By embedding NAG within an inertial
differential equation paradigm, they rigorously establish how varying the viscosity parameter alters
convergence rates and acceleration effects, bridging a crucial gap between continuous-time inertial
flow models and discrete-time iterative schemes. Expanding on this inertial dynamics perspective,
Wang and Peypouquet (2024) [441] focus specifically on strongly convex functions, where they de-
rive an exact convergence rate for NAG by constructing a novel Lyapunov function. Unlike previous
results that provided only upper-bound estimates for convergence, their approach offers a precise
characterization of NAG's asymptotic behavior, reinforcing its accelerated rate of $O(1/k^2)$ in smooth,
strongly convex settings. Their work strengthens the geometric interpretation of NAG as a dis-
cretization of a second-order differential equation with damping, further cementing its connection
to continuous-time optimization dynamics.

Despite the theoretical consensus on NAG's superiority in convex optimization, Hermant et. al. (2024) [442] present an unexpected empirical and theoretical challenge to this assumption. Their study systematically compares deterministic NAG with Stochastic Gradient Descent (SGD) under convex function interpolation, revealing cases where SGD exhibits superior practical performance despite lacking formal acceleration guarantees. Their findings raise fundamental questions about the practical advantages of momentum-based methods in data-driven scenarios, particularly when stochastic noise interacts with interpolation dynamics. Applying NAG beyond classical convex optimization, Alavala and Gorthi (2024) [443] integrate it into medical imaging reconstruction, specifically for Cone Beam Computed Tomography (CBCT). They develop a NAG-accelerated least squares solver (NAG-LS), demonstrating substantial improvements in computational efficiency and image reconstruction quality. Their results indicate that NAG's ability to mitigate error propagation in iterative reconstruction algorithms makes it particularly well-suited for inverse problems in medical imaging. From a generalization perspective, Li (2024) [444] formulates a unified momentum framework encompassing NAG, Polyak's Heavy Ball method, and other stochastic momentum algorithms. By introducing a generalized momentum differential equation, he rigorously dissects the trade-off between stability, acceleration, and variance control in momentum-based optimization. His framework provides a cohesive theoretical structure for understanding how momentum-based techniques interact with gradient noise, particularly in high-dimensional stochastic settings.

Beyond convexity, Gupta and Wojtowytsch (2024) [445] rigorously analyze NAG's performance
in non-convex optimization landscapes, a setting where standard acceleration techniques are often
assumed ineffective. Their research establishes conditions under which NAG retains acceleration
benefits even in the absence of strong convexity, highlighting how NAG's momentum interacts
with saddle points, sharp local minima, and benign non-convex structures. Their work provides
a crucial extension of NAG beyond convex functions, opening new avenues for its application in
deep learning and high-dimensional optimization. Gopal and Ghosh (2021) [653] presented a gen-
eral method for identifying natural numbers that satisfy reciprocity laws, encompassing quadratic,
cubic, quintic, and higher prime reciprocity cases. It extends classical reciprocity principles by
analyzing the structure of residues and nonresidues through a systematic transformation based on
modular differences. The study introduces the concept of degrees of freedom in residue groupings,
offering a new perspective on the patterns and relationships underlying higher-order reciprocity
laws. Meanwhile, Razzouki et. al. (2024) [446] compile a comprehensive survey of gradient-
based optimization methods, systematically comparing NAG, Adam, RMSprop, and other modern
optimizers. Their analysis delves into theoretical convergence guarantees, empirical performance
benchmarks, and practical tuning considerations, emphasizing how NAG's momentum-driven up-
dates compare against adaptive learning rate strategies. Their survey serves as an authoritative
reference for researchers seeking to navigate the landscape of momentum-based optimization al-
gorithms. Shifting towards hardware implementations, Wang et al. (2025) [447] apply NAG to
digital background calibration in Analog-to-Digital Converters (ADCs). Their study demonstrates
how NAG accelerates error correction algorithms in high-speed ADC architectures, particularly in mitigating nonlinear distortions and improving signal-to-noise ratios (SNRs). Their results provide compelling evidence that momentum-based optimization transcends software applications, finding practical utility in high-performance electronic circuit design.

To further explore empirical performance trade-offs, Naeem et. al. (2024) [448] conduct an exhaustive empirical evaluation of NAG, Adam, and Gradient Descent across various convex and non-convex loss functions. Their results highlight that while NAG accelerates convergence in many cases, it can induce oscillatory behavior in certain settings, necessitating adaptive momentum tuning to prevent divergence. Their findings offer practical insights into optimizer selection strategies, particularly in deep learning architectures where gradient curvature varies dynamically. Finally, Campos et. al. (2024) [449] extend NAG to optimization on Lie groups, a fundamental class of non-Euclidean geometries. By adapting momentum-based gradient descent methods to Lie algebra structures, they establish new convergence guarantees for optimization problems on curved manifolds, an area crucial to robotics, physics, and differential geometry applications. Their work signifies a major extension of NAG's applicability, proving its efficacy beyond Euclidean space.

### 17.3.2.2 Analysis of Nesterov Accelerated Gradient Descent (NAG)

Let $f: \mathbb{R}^d \to \mathbb{R}$ be a continuously differentiable function with a unique minimizer:

$$ \theta^* = \arg\min_{\theta \in \mathbb{R}^d} f(\theta). \qquad (17.54) $$

We assume the L-Lipschitz Continuity of the Gradient

$$ \|\nabla f(\theta) - \nabla f(\theta')\| \le L\|\theta - \theta'\|, \quad \forall \theta, \theta' \in \mathbb{R}^d. \qquad (17.55) $$

and the Strong Convexity of $f(\theta)$ with Parameter $m$

$$ f(\theta') \ge f(\theta) + \nabla f(\theta)^T (\theta' - \theta) + \frac{m}{2} \|\theta' - \theta\|^2, \quad \forall \theta, \theta' \in \mathbb{R}^d. \qquad (17.56) $$

The strong convexity assumption ensures that the Hessian satisfies:

$$ mI \le \nabla^2 f(\theta) \le LI, \quad \forall \theta \in \mathbb{R}^d. \qquad (17.57) $$

In Gradient Descent, Classical gradient descent updates follow:

$$ \theta_{t+1} = \theta_t - \eta \nabla f(\theta_t). \qquad (17.58) $$

This method achieves a linear convergence rate of $O(1/t)$ in the convex case. In Momentum-Based Gradient Descent, the momentum-based update rule is:

$$ v_t = \mu v_{t-1} - \eta \nabla f(\theta_t), \qquad (17.59) $$

$$ \theta_{t+1} = \theta_t + v_t. \qquad (17.60) $$

where $v_t$ is a velocity-like term accumulating past gradients. $\mu$ is the momentum coefficient. Momentum reduces oscillations and accelerates convergence but suffers from excessive oscillations in ill-conditioned problems. The Nesterov Accelerated Gradient (NAG) is A Look-Ahead Strategy. Instead of computing the gradient at $\theta_t$, NAG applies momentum first:

$$ \tilde{\theta}_t = \theta_t + \mu v_{t-1}. \qquad (17.61) $$

Then, the velocity update is performed using the gradient at $\tilde{\theta}_t$:

$$ v_t = \mu v_{t-1} - \eta \nabla f(\tilde{\theta}_t). \qquad (17.62) $$

Finally, the parameter update follows:

$$ \theta_{t+1} = \theta_t + v_t. \qquad (17.63) $$

The Interpretation of the Nesterov Accelerated Gradient (NAG) is
* **Look-Ahead Gradient Computation:** By computing $\nabla f(\tilde{\theta}_t)$ instead of $\nabla f(\theta_t)$, NAG effectively anticipates the next move, leading to improved convergence rates.

* **Adaptive Step Size:** The effective step size is modified dynamically, stabilizing the trajectory.

To find the Variational Formulation of NAG, We derive NAG from an auxiliary optimization problem that minimizes an upper bound on $f(\theta)$. Define a quadratic approximation at the look-ahead iterate $\tilde{\theta}_t$:

$$ \theta_{t+1} = \arg \min_{\theta} \left[ f(\tilde{\theta}_t) + \nabla f(\tilde{\theta}_t)^T (\theta - \tilde{\theta}_t) + \frac{1}{2\eta} \|\theta - \tilde{\theta}_t\|^2 \right]. \quad (17.64) $$

Solving for $\theta_{t+1}$:

$$ \theta_{t+1} = \tilde{\theta}_t - \eta \nabla f(\tilde{\theta}_t). \quad (17.65) $$

This derivation justifies why NAG achieves adaptive step-size behavior. We analyze the convergence properties and Optimality Rate under convexity assumptions of Gradient Descent (GD). For gradient descent:

$$ f(\theta_t) - f(\theta^*) = O\left(\frac{1}{t}\right). \quad (17.66) $$

This is suboptimal in large-scale settings. Regarding the NAG Convergence Rate, for strongly

<figure><img src="image_227.png" /><figcaption>Figure 17.3: Trajectory of Nesterov Accelerated Gradient (NAG) on $f(x, y) = x^2 + y^2$</figcaption></figure>

convex $f(\theta)$:

$$ f(\theta_t) - f(\theta^*) = O\left(\frac{1}{t^2}\right). \quad (17.67) $$
This improvement is due to the momentum-enhanced look-ahead updates. We need to do the Lyapunov Analysis for Stability. Define the Lyapunov function:

$$V_t = f(\theta_t) - f(\theta^*) + \frac{\gamma}{2} \|\theta_t - \theta^*\|^2. \qquad (17.68)$$

Here, $\gamma, \delta > 0$ are parameters chosen to ensure $V_t$ is non-increasing. We analyze $V_{t+1} - V_t$ to show it is non-positive. Expanding $V_{t+1}$:

$$V_{t+1} = f(\theta_{t+1}) - f(\theta^*) + \frac{\gamma}{2} \|\theta_{t+1} - \theta^*\|^2 + \frac{\delta}{2} \|v_{t+1}\|^2. \qquad (17.69)$$

Using strong convexity:

$$f(\theta_{t+1}) \le f(\theta_t) + \nabla f(\theta_t)^T (\theta_{t+1} - \theta_t) + \frac{L}{2} \|\theta_{t+1} - \theta_t\|^2. \qquad (17.70)$$

Since $\theta_{t+1} = \theta_t + v_t$, we substitute:

$$f(\theta_{t+1}) \le f(\theta_t) + \nabla f(\theta_t)^T v_t + \frac{L}{2} \|v_t\|^2. \qquad (17.71)$$

Now, using $v_t = \mu v_{t-1} - \eta \nabla f(\tilde{\theta}_t)$, we analyze the term $\|\theta_{t+1} - \theta^*\|^2$:

$$\|\theta_{t+1} - \theta^*\|^2 = \|\theta_t - \theta^* + v_t\|^2. \qquad (17.72)$$

Expanding:

$$\|\theta_{t+1} - \theta^*\|^2 = \|\theta_t - \theta^*\|^2 + 2(\theta_t - \theta^*)^T v_t + \|v_t\|^2. \qquad (17.73)$$

Similarly, we expand $\|v_{t+1}\|^2$:

$$\|v_{t+1}\|^2 = \|\mu v_t - \eta \nabla f(\tilde{\theta}_{t+1})\|^2. \qquad (17.74)$$

Expanding:

$$\|v_{t+1}\|^2 = \mu^2 \|v_t\|^2 - 2\mu\eta v_t^T \nabla f(\tilde{\theta}_{t+1}) + \eta^2 \|\nabla f(\tilde{\theta}_{t+1})\|^2. \qquad (17.75)$$

We have to choose $\gamma, \delta$ to Ensure Descent. To ensure $V_{t+1} \le V_t$, we require:

$$V_{t+1} - V_t \le 0. \qquad (17.76)$$

After substituting the above expansions and simplifying, we obtain a sufficient condition:

$$\gamma \ge \frac{L}{\eta}, \quad \delta \ge \frac{1}{\eta}. \qquad (17.77)$$

Choosing $\gamma, \delta$ appropriately, we conclude:

$$V_{t+1} \le V_t \qquad (17.78)$$

which proves the global stability of NAG. In conclusion, since $V_t$ is non-increasing and lower-bounded (by 0), it converges, which implies that $\theta_t \to \theta^*$ and the NAG iterates remain bounded. Hence, we have rigorously proven the global stability of Nesterov's Accelerated Gradient (NAG). For Practical Considerations, we need to have:

* **Choice of $\mu$**: Optimal momentum is $\mu = 1 - O(1/t)$.

* **Adaptive Learning Rate**: Choosing $\eta = O(1/L)$ ensures convergence.
### 17.3.2.3 Python Code to generate Figure 17.3

We'll again use the simple quadratic test function

$$f(x, y) = x^2 + y^2 \tag{17.79}$$

to visualize the optimization trajectory.

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

# Objective function
def f(x, y):
    return x**2 + y**2

# Gradient of the function
def grad_f(x, y):
    return np.array([2*x, 2*y])

# NAG parameters
lr = 0.1      # learning rate
mu = 0.9      # momentum
iterations = 20
x, y = 3.0, 3.0  # starting point

# Initialize velocity
v = np.array([0.0, 0.0])

# Store trajectory
trajectory = [(x, y)]

# NAG updates
for i in range(iterations):
    # Lookahead position
    lookahead = np.array([x, y]) - mu * v
    # Compute gradient at lookahead
    gradient = grad_f(lookahead[0], lookahead[1])
    # Update velocity
    v = mu * v + lr * gradient
    # Update position
    x, y = np.array([x, y]) - v
    trajectory.append((x, y))

trajectory = np.array(trajectory)

# Plotting
X = np.linspace(-4, 4, 100)
Y = np.linspace(-4, 4, 100)
X, Y = np.meshgrid(X, Y)
Z = f(X, Y)

plt.figure(figsize=(8,6))
plt.contour(X, Y, Z, levels=30, cmap='plasma')
plt.plot(trajectory[:,0], trajectory[:,1], 'ro-', label='NAG Path')
plt.title("Nesterov Accelerated Gradient (NAG) Trajectory on $f(x, y) = x^2 + y^2$")
```
48 plt.xlabel("x")

49 plt.ylabel("y")

50 plt.legend()

51 plt.show()

17.3.3 Adam (Adaptive Moment Estimation) Optimizer

17.3.3.1 Literature Review of Adam (Adaptive Moment Estimation) Optimizer

Kingma and Ba (2014) [176] introduced the Adam optimizer. It presents Adam as an adaptive
gradient-based optimization method that combines momentum and adaptive learning rate tech-
niques. The authors rigorously prove its advantages over traditional optimizers such as SGD and
RMSProp. Reddy et. al. (2019) [177] analyzed the convergence properties of Adam and identified
cases where it may fail to converge. The authors propose AMSGrad, an improved variant of Adam
that guarantees better theoretical convergence behavior. Jin et. al. (2024) [178] introduced MI-
Adam (Multiple Integral Adam), which modified Adam's update rules to enhance generalization.
The authors theoretically and empirically demonstrate its effectiveness in avoiding sharp minima.
Adly et. al. (2024) [179] proposed EXAdam, an improvement over Adam that uses cross-moments
in parameter updates. This leads to faster convergence while maintaining the adaptability of
Adam. Theoretical derivations show improved variance reduction in updates. Liu et. al. (2024)
[180] provided a rigorous mathematical proof of convergence for Adam when applied to linear in-
verse problems. The authors compare Adam's convergence rate with standard gradient descent and
prove its efficiency in noisy settings. Yang (2025) [181] generalized Adam by introducing a biased
stochastic optimization framework. The authors show that under specific conditions, Adam's bias
correction step is insufficient, leading to poor convergence on strongly convex functions. Park and
Lee (2024) [182] developed SMMF, a novel variant of Adam that factorizes momentum tensors, re-
ducing memory usage. Theoretical bounds show that SMMF preserves Adam's adaptability while
improving efficiency. Mahjoubi et al. (2025) [183] provided a comparative analysis of Adam, SGD,
and RMSProp in deep learning models. It demonstrates scenarios where Adam outperforms other
methods, particularly in high-dimensional optimization problems. Seini and Adam (2024) [184]
examined how Adam's optimization framework can be adapted to human-AI collaborative learning
models. The paper provides a theoretical foundation for integrating Adam into AI-driven education
platforms. Teessar (2024) [185] discussed Adam's application in survey and social science research,
where adaptive optimization is used to fine-tune questionnaire analysis models. This highlights
Adam's versatility outside deep learning.

17.3.3.2 Analysis of Adam (Adaptive Moment Estimation) Optimizer

The Adaptive Moment Estimation (Adam) optimizer can be considered a sophisticated, hybrid
optimization algorithm combining elements of momentum-based methods and adaptive learning
rate techniques, which is why it has become a cornerstone in the optimization of complex ma-
chine learning models, particularly those used in deep learning. Adam's formulation is centered
on computing and using both the first and second moments (i.e., the mean and the variance) of
the gradient with respect to the loss function at each parameter update. This process effectively
adapts the learning rate for each parameter, based on its respective gradient's statistical properties.
The moment-based adjustments provide robustness against issues such as poor conditioning of the
objective function and gradient noise, which are prevalent in large-scale optimization problems.

We aim to minimize a stochastic objective function $f(\theta)$, where $\theta \in \mathbb{R}^d$ represents the parameters of the model. The optimization problem is:

$$
\theta^* = \arg \min_{\theta} \mathbb{E}[f(\theta; \xi)] \qquad (17.80)
$$
where $\xi$ is a random variable representing the stochasticity (e.g., mini-batch sampling in deep learning). The Adam optimizer maintains that the first moment estimate (exponentially decaying average of gradients) is given by:

$$m_t = \beta_1 m_{t-1} + (1 - \beta_1) g_t \qquad (17.81)$$

where $g_t = \nabla_{\theta}f(\theta_{t-1}; \xi_t)$ is the stochastic gradient at time $t$, and $\beta_1 \in [0, 1)$ is the decay rate. The second moment estimate (exponentially decaying average of squared gradients) is given by:

$$v_t = \beta_2 v_{t-1} + (1 - \beta_2) g_t^2 \qquad (17.82)$$

where $\beta_2 \in [0, 1)$ is the decay rate, and $g_t^2$ denotes element-wise squaring. The bias-corrected estimates are:

$$\hat{m}_t = \frac{m_t}{1 - \beta_1^t}, \quad \hat{v}_t = \frac{v_t}{1 - \beta_2^t} \qquad (17.83)$$

The parameter update rule is:

$$\theta_t = \theta_{t-1} - \eta \frac{\hat{m}_t}{\sqrt{\hat{v}_t + \epsilon}} \qquad (17.84)$$

where $\eta$ is the learning rate, and $\epsilon > 0$ is a small constant for numerical stability. To rigorously analyze Adam, we impose the following assumptions. The gradient $\nabla_{\theta}f(\theta)$ is Lipschitz continuous with constant $L$:

$$\| \nabla_{\theta} f(\theta_1) - \nabla_{\theta} f(\theta_2) \| \le L \| \theta_1 - \theta_2 \| \qquad (17.85)$$

The stochastic gradients $g_t$ are bounded almost surely:

$$\|g_t\|_{\infty} \le G \qquad (17.86)$$

The second moments of the gradients are bounded:

$$\mathbb{E}[\|g_t\|^2] \le \sigma^2 \qquad (17.87)$$

The feasible region $\Theta$ is bounded with diameter $D$:

$$\|\theta_1 - \theta_2\| \le D, \quad \forall \theta_1, \theta_2 \in \Theta \qquad (17.88)$$

The decay rates $\beta_1$ and $\beta_2$ satisfy $0 \le \beta_1, \beta_2 < 1$, and $\beta_1 < \beta_2$. We analyze Adam in the online optimization framework, where the loss function $f_t(\theta)$ is revealed sequentially. The goal is to bound the regret:

$$R(T) = \sum_{t=1}^{T} f_t(\theta_t) - \min_{\theta} \sum_{t=1}^{T} f_t(\theta) \qquad (17.89)$$

The regret can be decomposed as:

$$R(T) = \underbrace{\sum_{t=1}^{T} f_t(\theta_t) - \sum_{t=1}^{T} f_t(\theta^*)}_{\text{Regret due to optimization}} + \underbrace{\sum_{t=1}^{T} f_t(\theta^*) - \min_{\theta} \sum_{t=1}^{T} f_t(\theta)}_{\text{Regret due to stochasticity}} \qquad (17.90)$$

Regarding the Boundedness of $\hat{m}_t$ and $\hat{v}_t$, using the boundedness of $g_t$, we can show:

$$\|\hat{m}_t\|_{\infty} \le \frac{G}{1 - \beta_1}, \quad \|\hat{v}_t\|_{\infty} \le \frac{G^2}{1 - \beta_2} \qquad (17.91)$$

The bias-corrected estimates satisfy:

$$\mathbb{E}[\hat{m}_t] = \mathbb{E}[g_t], \quad \mathbb{E}[\hat{v}_t] = \mathbb{E}[g_t^2] \qquad (17.92)$$