The Foundation of BioGA

Multi-Objective Genetic Algorithm for Genomic Data Optimization

Dany Mukesha

2026-10-17

Introduction

The BioGA utilizes multi-objective Genetic Algorithm (GA) specifically tailored for the optimization and analysis of genomic data. A detailed mathematical and theoretical analysis of the algorithms implemented within this package is described below; covering all major components of the evolutionary process with formal mathematical notation, proofs, and a discussion of key properties.

The Genetic Algorithm Framework

The BioGA package employs a standard generational Genetic Algorithm structure, modified for multi-objective optimization. This framework involves the iterative evolutionary cycle of initialization, evaluation, selection, reproduction (crossover and mutation), and replacement.

Mathematical Representation

Let us define the core elements of the search space and the algorithm’s state:

• $G = (V, E)$ is the gene network, where $V = \{g_1, g_2, \ldots, g_n\}$ are $n$ genes.
• $X \in \mathbb{R}^{n \times m}$ is the genomic data matrix, consisting of $n$ genes across $m$ samples.
• $P_t \in \mathbb{R}^{p \times n}$ is the population matrix at generation $t$, containing $p$ individuals (rows) and $n$ gene values (columns).

The evolution from one generation to the next, $P_t \to P_{t+1}$, is described by the following functional relationship:

$$P_{t+1} = R(M(C(S(P_t, f(P_t)), X)), P_t, f(P_t))$$

Where the functions represent the core operators:

• $f$: Fitness Evaluation function.
• $S$: Selection operator (NSGA-II inspired).
• $C$: Crossover operator (SBX).
• $M$: Mutation operator (Adaptive).
• $R$: Replacement operator (Elitism + Diversity Preservation).

Population Initialization

The initialization phase creates the starting population, $P_0$, by drawing random samples directly from the measured genomic data $X$.

Mathematical Formulation

Given the genomic data $X \in \mathbb{R}^{n \times m}$ and a population size $p$, the value of gene $j$ in the $i$-th individual of the initial population is set as:

$$P_0[i,j] = X[j, k] \quad \text{where} \quad k \sim \text{Uniform}\{1, \ldots, m\}$$ This means that for each element $(i, j)$ in $P_0$, a sample $k$ is chosen uniformly at random from the $m$ available samples, and the corresponding gene expression value $X[j, k]$ is used.

Initialization with Clustering:

If a clustering breakdown is provided (e.g., $C$ is the set of clusters), the initialization respects the structure: $$P_0[i,j] = X[j, k] \quad \forall j \in c, \quad k \sim \text{Uniform}\{1, \ldots, m\}$$ for all genes $j$ belonging to the same cluster $c$.

Statistical Properties

1. Distribution Preservation: This method ensures that the initial population individuals maintain the original statistical distribution of expression values for each gene $j$.
2. Structural Integrity: If clustering is used, the method preserves the pre-defined cluster structure within the initial population.
3. Expected Value: The expected value of any element $P_0[i,j]$ over the initialization process is the mean of the observed expression for gene $j$: $$\mathbb{E}[P_0[i,j]] = \mu_j = \frac{1}{m}\sum_{k=1}^m X[j,k]$$

Fitness Evaluation

The BioGA package uses a two-objective formulation. Both objectives are minimised. The package supports two representations, and the meaning of the objectives depends on which one is used.

Real representation

In the real representation an individual $i$ is a continuous expression profile $P_i \in \mathbb{R}^n$.

Objective 1: Data fidelity ($f_1$). The squared Euclidean distance between the individual and every observed sample. It is minimised when each gene value equals the per-gene mean $\mu_j$:

$$f_1(i) = \sum_{j=1}^n \sum_{k=1}^m (X_{jk} - P_{ij})^2$$

Objective 2: Parsimony ($f_2$). The mean absolute gene value (an $L_1$ penalty). It is minimised by driving gene values toward zero:

$$f_2(i) = \frac{1}{n}\sum_{j=1}^n |P_{ij}|$$

Binary representation

In the binary representation an individual $i$ is a selection indicator $P_i \in \{0,1\}^n$ and class labels $y \in \{1,\dots,c\}^m$ are required.

Objective 1: Predictive signal ($f_1$). For each gene $j$ an ANOVA $F$-statistic measures how strongly it separates the classes: $$F_j = \frac{\mathrm{SSB}_j/(c-1)}{\mathrm{SSW}_j/(m-c)}$$ where $\mathrm{SSB}_j$ and $\mathrm{SSW}_j$ are the between- and within-class sums of squares. The objective is the negated mean $F_j$ over the selected genes (so that better-separating subsets have a smaller value): $$f_1(i) = -\frac{1}{|S_i|}\sum_{j \in S_i} F_j, \qquad S_i = \{j : P_{ij} > 0.5\}$$ with $f_1(i) = 0$ when no gene is selected.

Objective 2: Sparsity ($f_2$). The fraction of selected genes: $$f_2(i) = \frac{|S_i|}{n}$$

Because $f_1$ rewards including informative genes while $f_2$ rewards including few genes, the two objectives conflict and the solution is a Pareto front rather than a single point.

Selection (NSGA-II)

The selection mechanism is based on the principles of the Non-dominated Sorting Genetic Algorithm II (NSGA-II), ensuring that evolutionary pressure drives the population toward the Pareto front.

Domination Criteria

Individual $i$ dominates individual $j$ (denoted $i \succ j$) if and only if:

1. Individual $i$ is no worse than $j$ on all objectives, and
2. Individual $i$ is strictly better than $j$ on at least one objective.

Formally: $$i \succ j \quad \text{iff} \quad \forall k \in \{1, 2\}: f_k(i) \leq f_k(j) \quad \text{and} \quad \exists k \in \{1, 2\}: f_k(i) < f_k(j)$$

Proof of Partial Order

The two-objective domination relation $(\succ)$ defines a partial order on the solution set because it satisfies the following properties:

1. Irreflexive: No individual dominates itself. If $i = j$, the second condition ($\exists k: f_k(i) < f_k(j)$) is false.
2. Antisymmetric Property: If $i \succ j$, then $j$ cannot dominate $i$. If $i$ is strictly better on at least one objective, $j$ must be strictly worse ($\forall k: f_k(j) \not\leq f_k(i)$).
3. Transitive: If $i \succ j$ and $j \succ k$, then $i \succ k$. This follows directly from the transitivity of the combined $\leq$ and $<$ relations.

Front Construction

The population is sorted into non-dominated fronts:

1. Domination Metrics: For every individual, we compute its domination count (number of individuals that dominate it) and its dominated set (set of individuals it dominates).
2. First Front ($\mathcal{F}_1$): This front consists of all individuals whose domination count is $0$. These are the current non-dominated solutions.
3. Subsequent Fronts: Individuals in $\mathcal{F}_t$ are temporarily removed. The domination counts of the individuals they dominated are reduced by one. The next front, $\mathcal{F}_{t+1}$, consists of all remaining individuals whose updated domination count is $0$. This process continues until the entire population is sorted.

Theorem on Computational Complexity

The non-dominated sorting and front construction algorithm terminates in quadratic time complexity relative to the population size.

Theorem: The front construction algorithm terminates in $O(p^2o)$ time, where $p$ is the population size and $o$ is the number of objectives.

Proof Sketch:

• The pairwise domination check between any two individuals is an $O(o)$ operation.
• The computation of the domination count and dominated set for all individuals requires $O(p^2o)$ time, as it involves $\binom{p}{2}$ comparisons.
• Once these metrics are computed, the iterative process of assigning individuals to fronts requires $O(p)$ time per front, leading to an overall complexity dominated by the $O(p^2o)$ initial comparison phase.

Crossover (Simulated Binary Crossover - SBX)

The Simulated Binary Crossover (SBX) is a prominent operator for real-coded GAs, designed to mimic the convergence properties of binary crossover while operating directly on continuous variables.

Given two parent individuals, $x, y \in \mathbb{R}^n$, the SBX operator generates two complementary offspring vectors. For each gene $j$:

With crossover probability $p_c$, the offspring genes are determined: 1. A random number $u$ is drawn: $u \sim \text{Uniform}(0,1)$. 2. The spreading factor $\beta$ is computed based on the distribution index $\eta$: $$\beta = \begin{cases} (2u)^{1/(\eta+1)} & \text{if } u \leq 0.5 \\ \left(\frac{1}{2(1-u)}\right)^{1/(\eta+1)} & \text{otherwise} \end{cases}$$ 3. The two complementary gene values are calculated: $$c_1 = 0.5[(1+\beta)x_j + (1-\beta)y_j], \qquad c_2 = 0.5[(1-\beta)x_j + (1+\beta)y_j]$$ and one of them is assigned to the offspring.

Otherwise (with probability $1-p_c$) the gene is copied unchanged from one of the parents. In the binary representation uniform crossover is used instead: each gene is copied from a randomly chosen parent.

Properties of SBX

1. Mean Preservation: The two complementary children satisfy $c_1 + c_2 = x_j + y_j$, so the expected value of an offspring gene (with the child chosen uniformly at random) is the parents’ average: $\mathbb{E}[z_j] = \frac{x_j + y_j}{2}$.
2. Distribution Control: The distribution index $\eta$ controls the proximity of the offspring to the parents.
3. Low $\eta$ (e.g., $\eta \to 0$): Produces a wide spread and approaches uniform crossover, generating many extreme values.
4. High $\eta$ (e.g., $\eta \to \infty$): Clusters offspring tightly around the parents, approaching no crossover.

Mutation

The mutation operation introduces random small changes to the offspring for exploration. The BioGA uses an adaptive mutation rate and incorporates network constraints.

For each gene $j$, with adaptive probability $p_m(t)$: $$p_m(t) = p_0\left(1 + 0.5\frac{t}{T}\right)$$ where $p_0$ is the initial rate, $t$ is the current generation, and $T$ is the total number of generations. This rate increases linearly over time to enhance exploration in later stages.

In the real representation the change applied to gene value $z_j$ is:

1. A Gaussian perturbation is generated: $\Delta_j \sim N(0, \sigma^2)$.

2. If a gene network $N$ is provided (adjacency matrix $N_{jk}$), the perturbation is damped by the gene’s relative connectivity. With the degree $d_j = \sum_{k \neq j}|N_{jk}|$ and $d_{\max} = \max_j d_j$: $$z_j \leftarrow z_j + \Delta_j\left(1 - \frac{d_j}{d_{\max}}\right)$$ so highly connected genes mutate less.

3. If no network is provided, standard Gaussian mutation is applied: $$z_j \leftarrow z_j + \Delta_j$$

In the binary representation mutation is a bit flip: with the adaptive probability the selected/not-selected state of a gene is inverted.

Properties of Adaptive Network Mutation

1. Adaptive Rate: The mutation rate increases with generation $t$, maintaining exploration later in the run.
2. Network Constraint Impact: The damping factor $(1 - d_j/d_{\max}) \in [0, 1]$ reduces the mutation magnitude for highly connected genes, preserving known network structure. It is a normalised, non-negative quantity, so mutation is never amplified.
3. Expected Change: Since $\mathbb{E}[\Delta_j] = 0$, the expected change remains zero: $\mathbb{E}[\Delta z_j] = 0$.
4. Variance Control: The variance of the change is $$\text{Var}(\Delta z_j) = \sigma^2\left(1 - \frac{d_j}{d_{\max}}\right)^2,$$ which shrinks monotonically with the gene’s connectivity.

Replacement

The replacement operator $R$ constructs the next population $P_{t+1}$ from the current population $P_t$ and the offspring $O_t$ using NSGA-II environmental selection.

1. Rank and crowding. The candidates are sorted into non-dominated fronts and a crowding distance is computed within each front (as in selection).
2. Elitist retention. The best $p - k_0$ parents, where $k_0$ is the replacement budget num_to_replace, are retained.
3. Offspring admission. The best $k_0$ offspring, ranked by the same (rank, crowding) criterion, fill the remaining slots.

Because $\text{rank}$ dominates $\text{crowding}$, non-dominated solutions are never discarded in favour of dominated ones; ties are broken towards less crowded (more diverse) solutions. The best solution with respect to the first objective is non-dominated by construction, so it is always retained when $k_0 < p$, giving monotonic non-increasing best $f_1$. If num_to_replace is zero the population is returned unchanged.

Convergence Analysis

BioGA is a heuristic. It does not provide a guarantee of finding the global Pareto front in a finite number of generations; the results below only outline the idealised conditions under which asymptotic convergence is expected, and should not be read as a practical guarantee.

Assumptions

For convergence results in MOEAs, the following are typically assumed:

1. Finite Search Space: While the gene values are continuous, they can be discretised in theory (or bounded in practice). The real-coded search space is not finite, so this assumption does not strictly hold for "real" mode.
2. Strictly Positive Mutation Probability: $p_m(t) > 0$ for all $t$, which enables ergodicity. This fails when mutation_rate = 0.
3. Elitism: The best solutions are preserved across generations. This holds only when num_to_replace < population_size.

Result (idealised)

Statement: Under the assumptions above, standard MOEA theory (e.g. Rudolph 1998, and analyses of NSGA-II) indicates that the population converges in probability to the true Pareto front as $t \to \infty$.

Sketch:

1. Preservation of optimality: Elitist replacement retains non-dominated solutions across generations.
2. Ergodicity: A strictly positive mutation probability lets any point be reached from any other in finitely many steps.
3. Convergence: Selection pressure toward non-dominated solutions, combined with ergodic variation, drives the population toward the Pareto front in the limit.

In practice, the useful convergence diagnostics are the per-generation best objective (history) and the final Pareto front; the theoretical guarantee above requires assumptions that real runs may violate.

Computational Complexity

Analyzing the complexity is crucial for understanding the algorithm’s performance on large genomic datasets. Key variables are defined:

• $p$ = population size
• $n$ = number of genes
• $m$ = number of samples
• $o$ = number of objectives ($o=2$ in BioGA)
• $T$ = number of generations

Component Complexities

GA Component	Complexity (Per Generation)	Dominating Variables
Initialization	$O(pn)$	$p, n$
Fitness Evaluation	$O(pn) \times O(m)$$\to$$O(pmn)$	$p, m, n$
Selection (NSGA-II sort)	$O(p^2o)$	$p$ (quadratic)
Crossover (SBX)	$O(pn)$	$p, n$
Mutation (Adaptive)	$O(pn)$	$p, n$
Replacement (environmental selection)	$O((p+q)^2o)$	$p, q$

Total Complexity

The total complexity for $T$ generations is dominated by the Fitness Evaluation and the Selection phases:

$$\text{Total Complexity} = O(T \cdot (pmn + p^2o))$$ Factoring $p$: $$\text{Total Complexity} = O(Tp(mn + po))$$ This result highlights that scalability is limited by the quadratic dependency on population size $p$ (due to sorting) and the linear dependency on the product of genes/samples $mn$ (due to fitness evaluation).

Mathematical Optimization Interpretation

From a purely mathematical perspective, the BioGA algorithm serves as a robust stochastic optimization solver for the following multi-objective problem:

$$\begin{aligned} \text{minimize } & (f_1(P), f_2(P)) \\ \text{subject to } & P \in \mathbb{R}^{p \times n} \end{aligned}$$

Core Characteristics that Favor GA:

1. Multi-Objective Nature: Explicitly seeking a set of trade-off solutions (the Pareto front).
2. High Dimensionality: The real search space is $p \times n$, making deterministic grid search intractable.
3. Non-Differentiable Objective: The sparsity objective $f_2$ in the binary representation is non-differentiable, so gradient-based methods do not apply.
4. Non-Convexity: The combined fitness landscape may be non-convex, so global search heuristics are useful to avoid poor local optima.

Biological Interpretation and Relevance

The mathematical operations within the BioGA framework have direct analogues in the biological context of genomic data analysis, ensuring the algorithm is biologically meaningful.

Mathematical Component	Biological Concept	Function
Population Initialization	Initial Biological Variability	Starting the search process by sampling observed expression levels.
Fitness ($f_1, f_2$)	Functional Efficacy & Parsimony	$f_1$ ensures the optimized solution is relevant to observed data; $f_2$ enforces biological simplicity (Minimal Regulatory Model).
Network Constraints	Gene-Gene Interactions	Incorporating known biological constraints (e.g., regulatory or metabolic pathways) to guide mutation.
Clustering	Co-expressed Gene Modules	Respecting known functional groupings during initialization.

Conclusion

The BioGA package implements a structured multi-objective evolutionary algorithm. Its components—NSGA-II rank and crowding for selection and environmental replacement, SBX (or uniform) crossover, and adaptive mutation with optional degree-based network damping—provide a coherent framework for trade-offs between data fidelity (or predictive signal) and model sparsity in genomic data optimisation.

Session Info

sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value
#>  version  R version 4.6.1 (2026-06-24)
#>  os       Ubuntu 24.04.5 LTS
#>  system   x86_64, linux-gnu
#>  ui       X11
#>  language en
#>  collate  C.UTF-8
#>  ctype    C.UTF-8
#>  tz       UTC
#>  date     2026-10-17
#>  pandoc   3.8.3 @ /opt/hostedtoolcache/pandoc/3.8.3/x64/ (via rmarkdown)
#>  quarto   NA
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package     * version date (UTC) lib source
#>  bookdown      0.48    2026-08-28 [1] RSPM
#>  bslib         0.12.0  2026-08-04 [1] RSPM
#>  cachem        1.1.0   2024-05-16 [1] RSPM
#>  cli           3.6.6   2026-04-09 [1] RSPM
#>  desc          1.4.3   2023-12-10 [1] RSPM
#>  digest        0.6.39  2025-11-19 [1] RSPM
#>  evaluate      1.0.5   2025-08-27 [1] RSPM
#>  fastmap       1.2.0   2024-05-15 [1] RSPM
#>  fs            2.1.0   2026-04-18 [1] RSPM
#>  htmltools     0.5.9   2025-12-04 [1] RSPM
#>  htmlwidgets   1.6.4   2023-12-06 [1] RSPM
#>  jquerylib     0.1.4   2021-04-26 [1] RSPM
#>  jsonlite      2.0.0   2025-03-27 [1] RSPM
#>  knitr         1.52    2026-09-06 [1] RSPM
#>  lifecycle     1.0.5   2026-01-08 [1] RSPM
#>  otel          0.2.0   2025-08-29 [1] RSPM
#>  pkgdown       2.2.1   2026-07-07 [1] RSPM
#>  R6            2.6.1   2025-02-15 [1] RSPM
#>  ragg          1.5.2   2026-03-23 [1] RSPM
#>  rlang         1.3.0   2026-07-05 [1] RSPM
#>  rmarkdown     2.32    2026-09-01 [1] RSPM
#>  sass          0.4.10  2025-04-11 [1] RSPM
#>  sessioninfo   1.2.4   2026-06-04 [1] RSPM
#>  systemfonts   1.3.2   2026-03-05 [1] RSPM
#>  textshaping   1.0.5   2026-03-06 [1] RSPM
#>  xfun          0.62    2026-10-16 [1] RSPM
#>  yaml          2.3.12  2025-12-10 [1] RSPM
#> 
#>  [1] /home/runner/work/_temp/Library
#>  [2] /opt/R/4.6.1/lib/R/site-library
#>  [3] /opt/R/4.6.1/lib/R/library
#> 
#> ──────────────────────────────────────────────────────────────────────────────