BSO was proposed by Shi [30] as a swarm-intelligence algorithm that models the divergent and convergent phases of human brainstorming rather than the foraging or social behaviors of animals. Every candidate solution is an idea, and every iteration partitions the population of ideas into clusters, generates new ideas from one or two clusters, and evaluates the resulting pool. The decisive structural difference with respect to classical evolutionary algorithms is that clustering carries much of the selection load: inferior individuals are not discarded outright but are retained as cluster members that keep seeding the generation, which preserves the diversity on which exploratory search depends [31], [32]. The original formulation performs clustering as a k-means-like partition around the best idea, and a selective strategy draws source ideas from better-performing clusters with higher probability.

In genetic-algorithm terms, the paradigm implements a soft, geography-preserving selection regime rather than a purely fitness-proportional one [33]. Its appeal for the present application is specific: Clusters of fusion configurations behave as coherent design lineages, and the inter-cluster generation step implements precisely the cross-family recombination that heterogeneous feature fusion requires, whereas a purely selection-driven search would collapse prematurely onto a single descriptor family. Comparative studies indicate that this organization is competitive with established evolutionary algorithms on multimodal and composition landscapes [31], which supports its adoption as the swarm backbone of a bilevel optimizer. The paradigm is therefore adopted unmodified at the level of its loop structure, with all modifications concentrated in the representation, the mutation operator, and the coordination between levels.

Feature extraction for EEG is organized around three complementary descriptor families, each instantiated by the case study of Section 6.2. Time-domain features summarize the morphology of the band-passed signal: Hjorth parameters (activity, mobility, complexity), channel variance, zero-crossing rate, and slope-based descriptors are inexpensive, tolerant of moderate artifacts, and historically foundational [4]. Frequency-domain features measure signal power in physiologically motivated bands; event-related desynchronization and synchronization over the sensorimotor cortex manifest as band-power modulations that are approximately log-normal, so logarithmic band powers serve as the standard representation [5]. Time– frequency features combine both views: CSP estimates spatial filters that maximize the ratio of between-class to within-class covariance [6], and filter-bank extensions repeat the estimation across overlapping sub-bands, so that subject-specific frequency foci emerge automatically [7]. When trials are scarce relative to dimension, regularized CSP and Riemannian covariance-based classifiers then markedly improve the stability of all three routes [8], [9].

Heterogeneous geometries are the outcome: A 40-dimensional time-domain block, a 32-dimensional frequency-domain block, and a 48-dimensional time–frequency block are estimated at different scales, exhibit strong within-block and weaker cross-block correlation, and contribute a total of 120 fused dimensions. Trial-wise concatenation without weighting leaves the downstream classifier exposed to scale mismatch, while hand-tuned weights ignore the interaction between fusion and regularization, an interaction resolved only when weighting and regularization are optimized in their correct hierarchical relationship. The fused representation is therefore treated as a designed artifact, optimized rather than fixed, which is the central methodological stance of this paper and what elevates the engineering task of Section 3.3 into the abstract bilevel framework of Section 3.1.

In the BCI setting, the empirical covariance of the fused feature space is almost always ill-conditioned, and the cause is the data regime itself: Trials per subject are of the same order as, or fewer than, the fused dimension, blocks are internally strongly correlated, and physiological artifacts inflate a few directions. An isotropically mutating optimizer therefore spends its evaluation budget on directions in which the objective is nearly flat, the deficiency that covariance matrix adaptation was designed to correct in evolution strategies [22], [26]. Such adaptation is most naturally organized on the cone of SPD matrices, whose geometry supports well-behaved means, distances, and tangent-space computations [34], [35] and which already underlies modern covariance-based BCI classifiers [9].

Two numerical safeguards are indispensable in measurement practice: 1) raw sample covariance matrices must be shrunk toward a well-conditioned target, for which the Ledoit–Wolf [36] estimator provides an analytically grounded rule, and 2) any update that leaves the SPD cone must be projected back, and eigenvalue clipping is the simplest such safeguard, with well-understood nearest-SPD interpretations [37]. HEF-BSO builds both safeguards directly into its mutation operator, so covariance awareness becomes a structural property of the search dynamics rather than an optional refinement. Under physiological noise, an additional benefit accrues: The follower's response mapping is evaluated with noise of physiological origin, and mutation distributions whose scale adapts to the local signal-to-noise geometry are markedly more robust under noisy evaluation [26], an issue to which the perturbation study of Section 7 returns quantitatively.

Bilevel optimization formalizes hierarchical decision problems in which a leader selects a decision x to minimize an upper-level objective F, while a follower, having observed x, responds by solving a nested lower-level problem in y [12], [13]. The optimistic formulation, adopted throughout this paper, is:

min ∗ F(x, y) s.t. y = y (x), (1)

where the follower's set of rational responses is given by

y (x) = {y ∈ (x): f(x, y) = miny ˜∈ (x) f(x, y ˜)}, (2)

and y*(x) denotes a selection from this set; pessimistic variants, which hedge against unfavorable selections, are not considered here. Implicit determination of the leader's feasible set through the follower's rationality is the source of the field's characteristic difficulties: Even linear–linear instances are strongly NP-hard [16], constraint qualifications fail systematically [14], and the follower's value function is generally discontinuous in x [15]. Multiobjective and set-valued extensions exist [18] but lie outside the present scope. The standard test suites, most prominently the SMD family [17], deliberately control these difficulties by pairing conflicting and non-conflicting leader–follower structures with convex, multimodal, or constrained lower levels and have become the de facto evaluation standard for evolutionary bilevel algorithms [19], [20]. Claims of progress on hierarchical problems must therefore rest on both final accuracy and anytime behavior on such suites, a requirement that shapes the experimental design of Section 6.

Whenever the lower-level problem is convex in y for every fixed x and suitable constraint qualifications hold, the follower's optimality can be replaced by its KKT conditions, which yields a single-level mathematical program with complementarity constraints [13], [38]. For lower-level inequality constraints g(x, y) ≤ 0 the stationarity condition reads

∇yf(x, y) + λ ∇yg(x, y) = 0, λ ≥ 0, (3)

together with primal feasibility, dual feasibility, and complementary slackness between the multipliers and the constraints; the full system is standard and is omitted for brevity. Direct attack on the reformulated problem by sequential quadratic programming defines the KKT-SQP route, which ranks among the strongest classical baselines when the lower level is well behaved [38]. Equally well documented are its weaknesses. At every feasible point, the complementarity constraints violate linear independence constraint qualifications, so SQP regularization schemes operate permanently on the boundary of admissibility. Moreover, nonconvex lower levels, as in the SMD instances with multimodal or quadratic followers, admit follower KKT points that need not be global optima, and the reformulation then solves a problem other than the one posed.

Population-based methods do not escape these issues for free, but degrade gracefully. When the follower's stationarity residual is treated as one term of a penalized leader objective, asymptotic correctness survives under mild conditions; an explicitly solved follower keeps validity under lower-level nonconvexity, at the price of an enlarged evaluation budget. HEF-BSO takes the second route, explicit follower solves coordinated by warm starts; the KKT-SQP baseline of Section 6.3 takes the first. Section 6.5 therefore probes this trade-off directly instead of presuming either reformulation philosophy superior.

The EEG measurement task instantiates the abstract program as follows. The leader variable x = (w, m) concatenates a vector of per-band fusion weights w in R³, one scalar gain for each of the time-domain, frequency-domain, and time–frequency blocks, with a soft selector mask m in [0, 1]¹² over twelve sub-band groups that attenuates or removes descriptors judged uninformative, giving p = 15 upper-level dimensions. Two regularization parameters per block compose the follower variable y in R⁶: A ridge penalty for the linear discriminant trained on the time- and frequency-domain descriptors and a CSP shrinkage intensity for the time–frequency block. For a given x, the lower-level problem minimizes the cross-validated training loss of the regularized classifier over y, and the upper-level objective is F(x, y) = 1 − κ, κ being Cohen's coefficient of agreement between predicted and actual imagery classes [39], so minimizing F maximizes chance-adjusted classification agreement.

Convexity of the lower level in y for fixed x holds in the practical operating range, because the regularized training loss is smooth and locally convex around the interior optimum; this justifies the warm-started inner solver of Section 4.5 without implying that the composite problem is easy. Piecewise smooth, mildly nonconvex, and noisy at the level of single evaluations, the composition of a resampling-based kappa estimate with a regularized inner fit is precisely the regime in which derivative-based leader updates are unreliable and in which population-based search with variance-scaled mutation is expected to pay off. The instance is thus representative of a broad class of hierarchical model-selection problems in which the leader designs a representation and the follower fits a model within it, and its structure, rather than any EEG-specific detail, is what the algorithmic machinery of the next section exploits.

In HEF-BSO, an evolving idea is a pair coupling a point in the space of fusion configurations with the local search geometry that surrounds it:

p d (4)

s = (w, Σ), w ∈ R , Σ ∈ S++.

The symbol w gathers the leader-side fusion parameters, namely the block weights and the selector mask of Section 3.3; the fused feature space has dimension d = 120; and an SPD matrix Σ parameterizes the mutation distribution of the idea. Weighted concatenation of the whitened blocks yields the fused feature vector handed to the classifier,

⊤ ⊤ ⊤ ⊤

z = [w z , w z , w z ] , (5)

1 time 2 freq 3 tf

where z_time, z_freq, and z_tf denote the per-block whitened descriptors, and the scalars w₁, w₂, w₃ are the corresponding block gains inherited from the leader variable. Carrying Σ inside the representation makes mutation covariance-aware rather than merely adaptive: The matrix is estimated from the dispersion of ideas within a cluster, projected onto the SPD cone, and then reused both as the mutation covariance of the cluster and as the cross-block coupling operator of Section 4.4. The redundancy with w, which already parameterizes the classifier input, is deliberate: It is what lets the search distribution encode block-scale and cross-block correlation information that a flat parameter vector cannot express.

The population of N = 50 ideas is partitioned at every iteration into k = 5 clusters in the w-space of fusion configurations. A k-means-like rule governs assignment [40]: Every idea is attributed to the cluster whose centroid is nearest,

ci = minj∈{1,…,k} ‖wi − μj‖, (6)

where c_i is read as the index attaining the minimum in (6), μ_j is the centroid, that is, the mean idea, of cluster j, and the norm is Euclidean in the normalized fusion space of Section 4.4. A single averaging pass, rather than iteration to convergence, updates the centroids, which keeps the clustering cost linear in population size and preserves the algorithm's anytime character; the deliberate under-optimization of the partition also leaves residual overlap between clusters, which the diversity analysis of [32] marks as beneficial for sustained exploration. Once assigned, each cluster accumulates the statistics consumed by the generation phase: its best idea, its centroid, and the empirical covariance of its members' w-vectors.

Ranking of clusters by their best evaluated objective precedes generation, and the selective strategy of the BSO paradigm biases the choice of source clusters toward better-ranked ones with a probability that decays with rank [30]. The clustering phase thereby furnishes both the topology of the search, that is, which ideas lie near which, and the geometry of mutation, that is, how ideas within a lineage disperse. This organization substitutes for much of the selection pressure that a canonical genetic algorithm would exert through tournament or proportional selection [33], while the population remains intact as a reservoir of partially evaluated design hypotheses.

Mutation of existing ideas with cluster-level Gaussian distributions produces new ideas. For an idea i belonging to cluster c_i, the mutation step is

wi = wi + η ⋅ (0, Σc ), (7)

where η is a global step-size schedule decaying from 0.5 to 0.05 over the run and the estimated covariance of cluster c_i serves as the mutation covariance. Under the reference configuration, a cluster holds roughly ten ideas, far fewer than the parameter dimension p, so the raw sample covariance is rank-deficient and violently ill-conditioned; the estimate is therefore shrunk toward its diagonal before use:

Σ ˆ = (1 − ρ)Σ ˆ + ρ ⋅ diag(Σ ˆ), (8)

sh

with shrinkage intensity ρ = 0.3 in the reference configuration, a fixed choice made in the spirit of the Ledoit– Wolf rule [36]. The eigenvalue spectrum of the mutation covariance is clipped after every update:

−4 3

clip (λ) = min(max(λ, 10 ), 10 ), (9)

which confines the condition number of every mutation covariance and admits an interpretation as a projection onto the set of SPD matrices with bounded spectrum, in the spirit of nearest-SPD correction [37]. Geometrically, the operator admits a reading on the SPD manifold: mutation proceeds as a tangent-space step at the cluster covariance, shrinkage and clipping play the roles of retraction safeguards [34], [35], [41], and the stochastic update inherits the stability properties established for Riemannian stochastic gradient schemes [42]. Inter-cluster ideas complete the generation phase: with probability 0.2, a new idea mutates a convex combination of two cluster centroids, implementing the cross-lineage recombination that the brainstorming metaphor prescribes [30]. What distinguishes HEF-BSO from both isotropic BSO [31] and single-distribution covariance adaptation [26] is the combination of a rank-aware selective strategy, full-covariance within-cluster mutation, and occasional inter-cluster generation: Several covariances are maintained simultaneously, one per design lineage, at a fraction of the cost of independent CMA-ES instances.

Before evaluation, the fusion operator acts on every newly generated idea and carries two responsibilities: Renormalizing the blocks and coupling them through the cluster covariance. Each block of the fused vector is rescaled by renormalization so that its median absolute deviation over the current population equals one. The step is performed in the whitened space, downstream of the per-block whitening filters fitted on training data, so that scale mismatch among variances, logarithmic band powers, and CSP log-variances cannot dominate the distance computations used by clustering and the classifier alike.

The block-partitioned covariance Σ = (Σ_ab), with a and b ranging over the time, frequency, and time– frequency blocks, is then applied by the coupling step. Diagonal blocks scale within-block mutation; off-diagonal blocks transfer perturbations across blocks: for instance, an idea that shifts the frequency-domain weighting toward high beta is nudged in the time-domain block along the historical co-movement of the two blocks within its cluster. Herein lies the precise sense in which mutation is covariance-aware, because the off-diagonal structure encodes cross-block correlation that isotropic operators destroy. Whenever reliable off-diagonal estimates are absent, as in early iterations or small clusters, the shrinkage of Eq. (8) interpolates automatically toward block-diagonal behavior, so the operator degrades gracefully rather than failing abruptly. Because the same partitioned matrix serves as the mutation covariance in Eq. (7) and as the coupling operator here, the operator carries no hyperparameters beyond the clipping bounds and the shrinkage intensity.

Scoring a candidate fusion configuration demands a lower-level regularization solve, and re-solving it from scratch for every idea would squander the evaluation budget. HEF-BSO therefore keeps an archive of evaluated leader–follower pairs with their objective values and warm-starts each follower solve from the response of the nearest archived leader:

(0) ′ ∗ ′

y (x ) = y (xnn(x )), (10)

where x_nn(x') denotes the archived leader nearest to x' in the normalized fusion space. The follower solver is a projected quasi-Newton inner loop run to a fixed iteration allowance; because lower-level optima move smoothly with the leader in the neighborhood of the archive, the warm start typically reduces the inner solve to a small correction. The evaluation budget is split 60/40 between the levels: Sixty percent of all function evaluations go to leader-side candidate evaluations, including the follower solves they trigger, and forty percent to lower-level refinement, so that the follower response never bottlenecks the leader search. When a new leader falls within a tolerance radius of an archived one, the archive reuses the stored response exactly, a cache that further amortizes follower cost. Warm-start coordination of this kind appears to be the component that prior bilevel metaheuristics leave unused, and the ablation of Section 6.7 quantifies its contribution at roughly an order of magnitude of median gap on the SMD suite.

Algorithm 1 assembles these components into the complete procedure, which is anytime, because the incumbent best configuration is available after every iteration, and deterministic once the random seed is fixed. Cluster covariances are refreshed every five iterations, with an exponential decay of 0.9 applied to the previous estimate in between; the arrangement keeps the amortized linear-algebra cost bounded as quantified in Section 5.2 and explains the measured wall-clock overhead of 21–29% over a standard BSO implementation on identical hardware. The follower consumes forty percent of evaluations, the archive grows by at most N entries per iteration, and all per-iteration data structures are of size O(Nd + d²), so the memory footprint is dominated by the covariance archive.

Algorithm 1. HEF-BSO.

1: Input: instance (F, f, Y), population N = 50, clusters k = 5, 2: budget split 60/40, covariance decay 0.9, clipping [1e-4, 1e3], 3: shrinkage rho = 0.3, step size eta: 0.5 to 0.05, eigen period 5 4: Initialize w_i for i = 1..N by Latin hypercube sampling over X 5: for each leader x_i: solve lower level, obtain y*(x_i), evaluate F 6: set Sigma_i to the identity; archive (x_i, y*(x_i), F_i) 7: while evaluation budget not exhausted do 8: assign each w_i to cluster c_i by nearest centroid, eq. (6) 9: update centroids mu_j; rank clusters by best F; mark elite 10: for each idea i in random order do 11: pick source cluster biased toward the elite (selective rule) 12: with prob. 0.2: form inter-cluster idea from two centroids 13: mutate: w_i' = w_i + eta * N(0, shrunken cluster cov.), (7)-(8) 14: apply block fusion operator: renormalize blocks, couple via Sigma 15: re-project Sigma' onto SPD cone by eigenvalue clipping, (9) 16: warm-start follower: y0 = y*(x_nn) of nearest archived leader 17: run lower-level solver within the 40% follower budget 18: evaluate F(x_i', y*(x_i')); update the archive 19: accept w_i' if it improves the incumbent of its cluster 20: every 5 iterations: refresh cluster covariances with decay 0.9 21: update global best configuration and best objective F_best 22: return best fusion configuration s = (w_best, Sigma_best), F_best

Two implementation details merit emphasis: 1) the whitened, block-normalized space hosts all distance computations entering in Eq. (6), the archive lookup of Eq. (10), and the renormalization statistics, so the heterogeneous scales of Section 2.2 never contaminate the geometry of the search, and 2) the covariance refresh schedule decouples the frequency of eigendecompositions from the frequency of mutations: Mutations consume the most recent clipped covariance at negligible cost, whereas the expensive spectral updates are amortized over five iterations, an arrangement whose complexity is made explicit in Eq. (12).

Three assumptions govern the analysis of the optimistic bilevel problem of Section 3.1.

Assumption 1. The search domain X × Y is compact, and F is continuous and attains its infimum F*. Assumption 2. For every x, the follower solver returns a response within ε_t of the lower-level optimum, with ε_t → 0 as the follower budget grows, so that evaluations are asymptotically consistent.

Assumption 3. Uniform non-degeneracy holds for every mutation covariance; once clipped by Eq. (9), its spectrum lies in [10⁻⁴, 10³], so the smallest eigenvalue is bounded below by a strictly positive constant. The optimality potential and its one-step drift are defined by

(t) (t) ∗ (t+1) (t)

V = F − F , E[V ∣ F ] ≤ V − δ ⋅ P(A ∣ F ), (11)

best t t t

where the σ-algebra of events summarizes the history of the algorithm up to iteration t, and A_t denotes the event that iteration t produces an improvement of at least δ in the incumbent. Elitist archiving keeps the potential non-increasing; under Assumption 3, the mutation distribution of every cluster assigns positive probability to every open ball around the incumbent, so on any non-absorbed state the probability of A_t is bounded below by a positive constant, and Eq. (11) exhibits the potential as a non-negative supermartingale with strictly decreasing drift outside the absorbed set. The clustering dynamics are absorbed in the limit: Cluster covariances contract around incumbents because their estimation inputs are the accepted ideas, while the decayed covariance refresh and the inter-cluster generation of Section 4.3 maintain non-vanishing exploration, leaving absorption at the global incumbent as the only closed recurrent class. The standard supermartingale argument then yields convergence of the incumbent objective to F* with probability one, in the lineage of convergence arguments for elitist evolutionary selectors [43], [44].

Deliberately structural rather than quantitative, the argument furnishes no convergence rates, and Assumption 2 encodes the exact sense in which the follower budget must scale for the guarantee to hold: if the follower is starved, the leader converges to the optimum of a perturbed problem rather than of the posed one. Section 9 revisits both restrictions as validity threats, and the empirical budget-split sensitivity reported in Section 7.2 measures how sharply practical performance responds to violations of Assumption 2.

Every iteration spends its dominant costs on assigning N ideas to k clusters, updating the k cluster covariances of size d × d, and performing the follower solves triggered by accepted ideas. With single-pass assignment, the clustering cost is linear in the population and the dimension, and the rank-based selective strategy adds only a logarithmic sorting term. Each cluster covariance demands an eigendecomposition at cost cubic in d, but these are refreshed every five iterations and decayed multiplicatively in between, so the amortized contribution is d³/5 for the aggregate of the k clusters when k is treated as a small constant:

citer = O(Nd + Nlog⁡k + d /5). (12)

With N = 50, k = 5, and d = 120 under the reference configuration, the amortized eigen cost dominates at roughly 3.5 × 10⁵ floating-point operations per iteration [45], two orders of magnitude below the cost of the follower evaluations on the EEG instance. Memory is O(Nd + d²) for the population and the covariance archive. The wall-clock overhead relative to a standard BSO with identical population and budget, measured at 21–29% in Section 6, is attributable almost entirely to the amortized eigendecompositions and the archive bookkeeping; the follower budget, at forty percent of evaluations, is a reallocation rather than an addition, since any nested baseline must likewise evaluate the follower to score a leader. The accounting clarifies why covariance awareness comes cheap here: Expensive linear algebra is amortized across five iterations and shared across the k clusters, whereas a bank of independent covariance-adaptive instances would pay the cubic cost per distribution and per iteration [26].

The SMD family [17] is the standard controlled testbed for bilevel optimizers: Every instance is constructed with a known global optimum, and the family spans the structural factors that make bilevel search hard, namely inter-level conflict, lower-level multimodality, constraints, and nonconvexity. The study evaluates ten instances, SMD1 through SMD8, SMD10, and SMD12, whose properties Table 1 summarizes. SMD1 through SMD6 pair conflicting and non-conflicting leader–follower structures with convex, multimodal, or constrained lower levels; SMD7 and SMD8 add nonconvex follower geometries at larger dimensions; SMD10 couples a quadratic follower to a constrained leader; and SMD12 combines conflict, constraints, and lower-level nonconvexity at the largest scale. Coupling and variance settings were freshly drawn for this program: the interlevel coupling matrices carry instance-specific scale factors ranging from 0.4 to 3.0, and follower curvature spreads span an order of magnitude across the suite, so that both loose and tight coupling regimes are represented. The suite deliberately complements the EEG instance: SMD problems are smooth and exactly specified, isolating the optimizer's search dynamics from the estimation noise of physiological data.

*p and n denote upper- and lower-level dimensions; the coupling scale multiplies the entries of the interlevel coupling matrices, and the follower variance controls the curvature spread of the lower-level quadratic terms. All instances admit a known optimum with F* = 0.

Reporting relies on the optimality gap F(x_best) − F*, computed at the best feasible leader found, with feasibility with respect to the follower's optimality certified by the lower-level solve, following the evaluation practice of [19]. Performance appears per instance in Table 4 and in aggregate through the suite-wide median, with median convergence profiles in Fig. 1. Every algorithm draws on the same follower-certification machinery, so observed differences reflect leader-side search behavior rather than bookkeeping conventions.

A motor-imagery-style benchmark protocol governs the case study. Eleven volunteer subjects performed cue-guided binary, left- versus right-hand, imagery while sixteen EEG channels over the sensorimotor strip were recorded, band-pass filtered, and epoched into trials. Three feature blocks were extracted per trial, as summarized in Table 2: Forty time-domain descriptors (Hjorth activity, mobility, and complexity, channel variance, zero-crossing rate, and slope descriptors over eight channel pairs [4]); thirty-two frequency-domain descriptors (logarithmic band powers in the theta, mu/alpha, beta, and low-gamma bands over the same eight pairs [5]); and forty-eight time–frequency descriptors (filter-bank CSP log-variances over six sub-bands [6], [7]), each block whitened on the training subjects before weighted concatenation through Eq. 5. Validation follows leave-one-subject-out cross-validation: For each of the eleven folds, the bilevel pipeline, its follower classifier, and every decoder were fitted on the ten training subjects and scored once on the held-out subject, so the reported metrics measure transfer to an unseen, calibration-free subject rather than within-session discrimination. Balanced accuracy, Cohen's Kappa [39], and area under the receiver-operating curve constitute the reported metrics.

blocks via Eq. 5 and selector mask m *The fusion space dimension d = 120 is the dimension used in Sections 4 and 5.

The d of Sections 4 and 5 is the resulting 120-dimensional fusion space, and the bilevel instance is the one formalized in Section 3.3: The leader optimizes the three block gains and the twelve selector-mask entries on the training subjects, while the follower optimizes the six per-block regularization parameters under an inner resampling of those subjects. Because the kappa scale is chance-adjusted, it is more conservative than raw agreement; a leave-one-subject-out kappa of 0.51 corresponds to agreement well above chance for binary imagery classes and sits in the transfer-level range reported for motor-imagery pipelines in the BCI literature [1], [3]. Between-subject variability dominates BCI performance, so the per-subject results of Section 6.6 and the subject-level pooling of the hierarchical model carry the interpretive load.

The bilevel leg compares eleven configurations: The proposed optimizer, nine external rivals, and one ablated twin. KKT-SQP solves the single-level reformulation of Section 3.2 with a regularized trust-region SQP method [38]; bilevel PSO and bilevel DE nest the canonical swarm and differential operators inside an explicit follower solve [23], [25]; standard BSO applies the original brainstorming dynamics to the KKT reformulation [30]; and CMA-ES on KKT equips the same reformulation with covariance matrix adaptation [26]. Random search samples leader–follower pairs uniformly over the joint box, while nested CMA-ES runs covariance adaptation at both levels. BLEAQ, the bilevel evolutionary algorithm with quadratic approximations, evolves a leader population while refitting quadratic models of the lower-level optimal solution mapping [19], and a surrogate-assisted bilevel GA replaces exact follower solves with a radial-basis-function surrogate of the response mapping [19], [46]. HEF-BSO without covariance-aware mutation is the ablated variant whose mutation is isotropic under a matched step-size schedule, and on the decoding leg the fused pipeline meets EEGNet, a compact convolutional network built from depth-wise and separable convolutions [10], together with the ShallowConvNet and DeepConvNet temporal-convolution decoders of the same lineage [11], each retrained on the ten training subjects of every leave-one-subject-out fold. All bilevel methods share the evaluation budget, the default 60/40 split wherever nesting applies, and the follower solver of Section 4.5 wherever an explicit follower is required; Table 3 records the remaining settings, parameter choices follow the source publications wherever available, and no algorithm was tuned on the reported instances.

Two convolutional blocks; Trained per LOSO fold on dropout 0.25; early stopping ten subjects on fold validation Single temporal-conv block Trained per LOSO fold on with spatial filter; early ten subjects stopping Four convolutional blocks; Trained per LOSO fold on dropout 0.5; early stopping ten subjects *All bilevel configurations receive the same total evaluation budget; the 60/40 leader-follower split applies wherever an explicit or approximate follower is solved. Decoders share the trial-wise preprocessing chain of Section 6.2 and the same

Selection Of The Bilevel Comparison Set Follows The Architectural Questions That HEF-BSO Raises, And Each Entry Prices One Design Decision. KKT-SQP [38] Anchors The Derivative-Based Wing: It Prices The Strongest Reformulation Route, Converges Rapidly On Well-Conditioned Convex Followers, And Stalls On The Complementarity Degeneracy. Bilevel PSO, From The Nested Particle-Swarm Family [23], Supplies The Reference Nested Swarm, Whose Inertia-Weighted Velocities Yield Steady Early Progress But A Slow Terminal Phase. Bilevel DE Adds A Second Population-Based Operator With A Different Exploration Geometry [25], Effective Across Flat Regions Yet Weaker In Terminal Precision, And Standard BSO [30] Isolates The Brainstorming Backbone Itself, With Clustering-Driven Diversity As Its Strength And Isotropic

BLEAQ and the surrogate-assisted bilevel GA price the approximate-mapping school from two ends. BLEAQ [19] refits local quadratic models of the follower's optimal-response mapping around evaluated leaders, which collapses the cost of lower-level optimization on smooth, well-conditioned followers; its limitation is curvature fidelity, because a quadratic model loses accuracy exactly where the response bends under tight coupling or constraints. The radial-basis-function surrogate GA represents the coarser end of the same school [19], [46], trading approximation bias for even cheaper leader-side evaluations. Nested CMA-ES prices the strongest explicit two-level search, running covariance adaptation at both levels under the same budget split, but a single covariance per level, estimated across the whole population, cannot resolve the cluster-level, block-partitioned correlation structure of heterogeneous feature blocks. CMA-ES on KKT prices what one adapted covariance achieves when the follower is absorbed into the reformulation [26], random search supplies the floor of the comparison, certifying that the benchmark rewards structure rather than sampling volume, and the ablated twin serves as the direct counterfactual for H1, identical in representation, clustering, and warm-start coordination.

The decoding leg was selected around the transfer question that RQ3 poses. EEGNet [10] is the strongest compact convolutional decoder on small-sample EEG and regularizes aggressively through depth-wise and separable convolutions, which makes it the natural rival for a fusion pipeline that must survive subject-level transfer. ShallowConvNet and DeepConvNet [11] bracket the depth axis of the same temporal-convolution family: the shallow variant carries few parameters and stays robust on small cohorts, whereas the deep variant fits richer features and is expected to overfit part of an eleven-subject cohort. Together the fourteen configurations span reformulation against nesting, exact against approximate follower responses, single against cluster-level covariances, and handcrafted fusion against end-to-end decoding, so each architectural choice is priced by at least one competitor.

Each bilevel algorithm was executed 30 times per instance with independent seeds and a common budget of 550 leader iterations, corresponding to 25,000 leader-side evaluations at population 50, with the follower consuming forty percent of the total evaluation budget under the reference split wherever nesting applies. The per-instance test is the Pratt-adjusted Wilcoxon signed-rank test on paired final gaps across the 30 runs [47], which retains zero differences by the Pratt convention so that ties in terminal precision are not silently discarded; Holm's step-down procedure controls family-wise error at α = 0.05 across the ten instances within each rival comparison [48], applied within the recommendation framework for evolutionary-algorithm benchmarking [49], [50]. Probability-of-superiority effect sizes accompany every test, computed as the posterior probability that a randomly drawn HEF-BSO gap falls below the rival's and following the ordinal-analysis tradition of dominance statistics [51].

The decoding leg carries the primary analysis, a hierarchical Bayesian model fitted to the leave-one-subject-out balanced accuracies. Subject random intercepts and method fixed effects compose the mean structure, so each subject contributes one repeated-measures block measured under all four pipelines, and partial pooling shares strength across the eleven subjects instead of estimating each in isolation [52], [53]. Posterior means are reported with 95% highest-density intervals, and the decision quantities are the pairwise posterior probabilities P(HEF-BSO > rival | data) for the three decoders. Superiority is declared at a posterior probability above 0.95, a threshold calibrated against the conventional evidence scale [54], with values between 0.80 and 0.95 reported as credible but not decisive. The model was sampled with four Hamiltonian chains of 2,000 warmup and 4,000 retained draws each; convergence is certified by R-hat diagnostics, all below 1.01, and by posterior predictive checks on subject-level dispersion, in which simulated subject spreads match the observed ones within twelve percent. Fallbacks are fixed in advance: any R-hat at or above 1.01 triggers doubled warmup and resampling, and a failed dispersion check demotes the primary claim to subject-level sign tests with Holm correction.

The statistical structure mirrors the measurement design. Because every subject is measured under every method, the repeated-measures-within-subject pairing is preserved by the random-intercept structure rather than broken into independent two-sample comparisons; because eleven subjects support transfer claims only weakly under no pooling, the partial-pooling compromise shrinks unstable subject effects toward the cohort mean while letting strongly evidenced ones stand; and because instance-level gaps span five orders of magnitude, rank-based ordinal tests carry the bilevel leg where means and variances would be dominated by a single instance. Post-hoc analyses, namely the ablation, the coupling stress, and the perturbation study, are labeled exploratory and reported with medians and interquartile ranges rather than confirmatory tests. Reproducibility rests on fixed conventions: One random seed per run, Python 3.11 with NumPy and SciPy, the EEG preprocessing tooling of Section 6.2 on one workstation of sixteen cores and 64 GB of memory, and an acquisition and preprocessing chain, montage of sixteen channels, band-pass filtering, epoching, and per-block whitening, held fixed throughout. The measured wall-clock overhead of the covariance machinery ran 21 to 29% over a standard BSO implementation on identical hardware, and every configuration, seed list, and result table is archived with the code.

Table 4 lists per-instance medians and interquartile ranges of the final gap, and Fig. 1 shows the median convergence profiles across the suite. HEF-BSO attains the smallest median gap on seven of the ten instances, with a suite-wide median of 4.4 × 10⁻⁶ and an interquartile range across instances of 3.0 × 10⁻⁶ to 5.8 × 10⁻⁶. The strongest rival, nested CMA-ES, closes at a suite-wide median of 2.4 × 10⁻⁴ and wins SMD4 and SMD10, the two well-conditioned multimodal instances on which a single per-level covariance can concentrate the whole population in one exploitation basin. BLEAQ wins SMD1, whose loose coupling and convex follower flatten its quadratic response model, but degrades on the tightly coupled constrained instances, closing at 2.8 × 10⁻³ on SMD5, 4.1 × 10⁻³ on SMD6, and 9.6 × 10⁻³ on SMD12, a pattern that confines approximation-based search to loosely coupled territory. The surrogate-assisted bilevel GA sits mid-pack at a suite-wide median of 9.6 × 10⁻⁴, competitive on the unconstrained instances and behind the exact nested solvers elsewhere, while bilevel PSO, bilevel DE, standard BSO, CMA-ES on KKT, KKT-SQP, and random search trail in that order, with KKT-SQP stalling at 1.2 or higher on every instance.

rival, reached without early acceleration.

BLEAQ Bilevel PSO KKT-SQP 4.1e-7 [2.5e-7, 6.0e-7] 1.1e-3 [6.9e-4, 1.8e-3] 1.2 [0.7, 2.5] 6.3e-4 [3.5e-4, 1.0e-3] 1.9e-3 [1.3e-3, 3.1e-3] 1.6 [0.9, 3.3] 8.1e-4 [5.7e-4, 1.5e-3] 2.4e-3 [1.7e-3, 3.5e-3] 2.1 [1.2, 4.3] 9.4e-4 [5.2e-4, 1.4e-3] 3.1e-3 [2.2e-3, 5.4e-3] 2.4 [1.2, 4.3] 2.8e-3 [1.9e-3, 4.0e-3] 2.6e-3 [1.6e-3, 4.8e-3] 2.2 [1.3, 4.0] 4.1e-3 [2.4e-3, 7.5e-3] 2.9e-3 [1.8e-3, 4.1e-3] 2.6 [1.4, 4.9] 7.2e-4 [4.7e-4, 1.2e-3] 2.2e-3 [1.3e-3, 3.7e-3] 1.8 [0.9, 3.6] 6.8e-4 [4.6e-4, 1.1e-3] 2.8e-3 [1.9e-3, 4.3e-3] 2.9 [1.7, 5.6] 8.9e-4 [5.7e-4, 1.3e-3] 3.4e-3 [2.0e-3, 4.9e-3] 3.3 [1.7, 6.4] 9.6e-3 [6.8e-3, 1.7e-2] 4.7e-3 [3.4e-3, 7.6e-3] 3.8 [2.2, 6.6]

*Suite-wide medians: HEF-BSO 4.4e-6; nested CMA-ES 2.4e-4; HEF-BSO w/o cov. mutation 4.4e-4; BLEAQ 8.5e-4; surrogate-assisted bilevel GA 9.6e-4; bilevel PSO 2.7e-3; bilevel DE 5.9e-3; standard BSO 1.3e-2; CMA-ES on KKT 2.8e-2; KKT-SQP 2.3; random search 34.6.

The signed-rank battery consolidates the margins. Against nested CMA-ES, Holm-adjusted p-values fall between 1.4 × 10⁻⁷ and 0.038 on the eight instances where the HEF-BSO median is smaller, with SMD1 significant at 0.0031 despite BLEAQ's instance win, and the two reversals on SMD4 and SMD10 remain non-significant after correction at 0.062 and 0.11. Against BLEAQ, the adjusted p-values span 2.1 × 10⁻⁶ to 0.044 on the nine instances where the HEF-BSO median is smaller, while the SMD1 reversal is itself decisive in BLEAQ's favor at 0.008; against the remaining seven rivals, every instance-wise comparison rejects at 0.049 or below after correction. Probability-of-superiority values, derived from the pooled posterior over paired gaps, read 0.86 against nested CMA-ES, 0.88 against BLEAQ, 0.87 against the ablated twin, 0.90 against the surrogate GA, and 0.93 or above against every remaining rival, so the close-competitor margins stay honestly short of saturation. Anytime behavior splits the difference: The median curve of HEF-BSO crosses a gap of 10⁻² after roughly 380 leader iterations, whereas nested CMA-ES, BLEAQ, and bilevel PSO need 510 to 525, because the early-phase gaps that dominate a budget-integrated view are similar across population-based methods. The honest summary is a terminal-precision advantage of a factor of roughly fifty over the strongest HEF-BSO reaches a median final gap of 4.4 × 10⁻⁶; nested CMA-ES, BLEAQ, bilevel PSO, and KKT-SQP terminate at 2.4 × 10⁻⁴, 8.5 × 10⁻⁴, 2.7 × 10⁻³, and 2.3.

The primary outcome is the hierarchical Bayesian analysis of leave-one-subject-out balanced accuracy, summarized in Table 5. The fused pipeline attains a posterior mean of 75.5% with a 95% highest-density interval of [72.5, 78.5], against 74.1 [71.5, 76.7] for EEGNet, 71.1 [68.1, 74.1] for ShallowConvNet, and 68.8 [64.2, 73.4] for DeepConvNet. The posterior probability that HEF-BSO exceeds its rival is 0.90 against EEGNet, above the 0.80 credibility floor but short of the 0.95 decisiveness threshold, and 0.996 and 0.999 against ShallowConvNet and DeepConvNet. On the companion metrics, the ordering is preserved, with posterior mean kappa of 0.508 against 0.480, 0.419, and 0.377 for the decoders and posterior mean area under the curve of 0.824 against 0.802, 0.778, and 0.752. R-hat statistics stayed between 1.002 and 1.006 on every monitored quantity, and the posterior predictive check on subject-level dispersion passed, with simulated spreads within eleven percent of the observed ones.

*Subject random intercepts with method fixed effects; four Hamiltonian chains with 2,000 warmup and 4,000 retained draws each. HDI denotes the highest-density interval; the cohort-level posterior mean kappa and AUC accompany the balanced-accuracy scale.

Per-subject estimates in Table 6 and Fig. 2 expose the structure that the pooled means conceal. The fused pipeline leads on nine of eleven subjects, loses S3 and S9 to EEGNet by 1.7 and 2.1 points, and never falls below 67.5%, whereas DeepConvNet spans 55.6 to 80.1%, the largest subject-to-subject dispersion in the comparison, with its weakest transfers on the two lowest-trial-count subjects. Under partial pooling, the raw subject effects shrink toward the cohort mean by twenty-two percent on average: S5's raw +7.5-point deviation contracts to +5.9, S8's -8.2 to -6.4, and S3's +4.7 to +3.7. Two subjects, S4 and S8, remain more than five points below the cohort mean after shrinkage and are classified as non-generalizing, their fusion configurations failing to transfer even though the optimizer fits them well on the training side, while S5 and S7 remain above four points and carry the transfer claim. The shrinkage pattern indicates that the cohort, not any single subject, supports the pooled superiority statement, which is precisely the claim structure an eleven-subject transfer study can defend.

*The cohort row reports the mean across subjects with the between-subject standard deviation; cohort-level posterior means for Cohen's Kappa and AUC appear in Table 5.

HEF-BSO leads on nine of eleven subjects; EEGNet wins S3 and S9, and DeepConvNet shows the largest subject-to-subject dispersion.

Learned configurations retain their interpretability under transfer. On the training folds, the selector mask converged to sparse patterns, typically retaining the mu/alpha and beta groups of the frequency-domain block and two to three filter-bank sub-bands of the time–frequency block while attenuating the theta group, a pattern consistent with the physiological locus of motor-imagery discrimination [5], [7]. The learned block gains favored the time–frequency block most heavily on nine of eleven held-out subjects, in line with the comparative standing of filter-bank CSP in the BCI literature [3], [7]. The decoder side tells the complementary story: EEGNet's two subject wins come from clean, high-trial-count recordings where end-to-end feature learning has room to operate, and DeepConvNet's dispersion is the classical overfitting signature of a four-block convolutional stack trained on ten subjects. The fused pipeline wins by transferring a regularized, physiologically grounded representation rather than by re-learning one per cohort, which is the operational difference the leave-one-subject-out protocol was designed to expose.

The ablation strips each of the three principal components in isolation while budget, seeds, and all remaining parameters are held fixed; Table 7 summarizes the outcomes. Stripping covariance-aware mutation hurts most: the suite-wide median gap degrades from 4.4 × 10⁻⁶ to 4.4 × 10⁻⁴, a factor of 100, confirming that isotropic mutation squanders the budget in flat and correlated directions. Removing warm-start bilevel coordination degrades the median to 5.2 × 10⁻⁵, a factor of 11.8, with the effect concentrated on the larger-follower instances SMD8, SMD10, and SMD12, where cold-started inner solves waste the follower allowance, and dropping heterogeneous block weighting degrades the median to 3.6 × 10⁻⁵, a factor of 8.2, because without per-block renormalization the largest-variance block dominates the distance computations and the classifier input alike.

*Factors are rounded to one decimal place; the factor of 100 for the isotropic variant is exact to the reported precision.

The drift of the learned covariance spectra supplies the mechanism behind these factors. Across the suite, the median condition number of the cluster covariances falls from 2.6 × 10³ after the first refresh to 6.1 × 10² at the midpoint and 1.4 × 10² at termination, a monotone contraction that tracks the exploitation phase. On the EEG instance, the log-spread of block-level eigenvalues contracts from 2.9, 2.4, and 3.3 decades for the time, frequency, and time–frequency blocks respectively, at initialization to 0.8, 1.1, and 0.9 decades at convergence, while the cross-block share of spectral mass rises from 11 to 34%. The two trajectories encode complementary information: within-block contraction measures how sharply each design lineage has focused, and cross-block growth measures how much coupling the accepted ideas have accumulated. The isotropic ablation can express neither quantity, which is why its damage dwarfs the other two, and the block-weighting ablation destroys the scale-free basis on which the spread is computed, which explains its smaller but systematic effect.

Stable across the suite is the ordering of the three effects, covariance awareness first, warm-start coordination second, and block weighting third, and the ablated variants remain competitive with every external rival except nested CMA-ES and BLEAQ on their winning instances. On this suite, the component contributions are individually necessary and jointly sub-additive: no single component reproduces the terminal precision of the full algorithm, yet the damage inflicted by removing any one of them stays below the sum of the pairwise gaps. That sub-additivity is the expected signature of components that partially compensate through the shared covariance and archive machinery.

Interlevel coupling is the load-bearing wall of a bilevel solver, and deployed hierarchical problems rarely afford a well-fed follower, a concern documented since the earliest bilevel bibliographies [55], [56]. This subsection therefore stress-tests the coupling along two axes: the follower's share of the evaluation budget, swept over 30, 40, 50, and 60%, that is, follower-leader splits of 30/70, 40/60, 50/50, and 60/40, and the constraint-coupling matrices of the suite, scaled by factors from 0.25 to 4. Two quantities are tracked over 30 independent runs: the bilevel infeasibility rate, the fraction of final leader solutions whose returned follower response violates lower-level optimality beyond a tolerance, and the median leader-side gap. Nested CMA-ES, the strongest rival of Section 6.5, serves as a comparator under the identical protocol, and the design mirrors deployment, where follower evaluations form the expensive half of the measurement chain.

The protocol remains deliberately lean, attaching to any nested bilevel run. For each budget share, each of the ten SMD instances, and each of the 30 independent runs, HEF-BSO and nested CMA-ES are executed with the follower allowance fixed at that share, leaving the total budget of 25,000 leader-side evaluations untouched. A high-accuracy reference solve certifies lower-level optimality, and a final leader solution counts as bilevel infeasible when its follower response misses that reference beyond the tolerance. The constraint-coupling scale is then swept on a coarser grid, with the median leader-side gap recorded per value and BLEAQ added as a third arm because its quadratic response models are the theory's predicted failure point under tightened coupling.

Table 8 reports the budget sweep, and Fig. 3 plots the median leader-side gap against the follower share. At the reference 40% share, the row reproduces the suite-wide medians of Section 6.5, namely 4.4 × 10⁻⁶ for HEF-BSO and 2.4 × 10⁻⁴ for nested CMA-ES, anchoring the stress test to the main results. Starving the follower to thirty percent leaves HEF-BSO at 7.8% infeasibility with a median gap of 4.6 × 10⁻⁵, where nested CMA-ES reaches 24.6% and 8.4 × 10⁻⁴, while feeding it to sixty percent buys little further accuracy, 5.1 × 10⁻⁶ against 2.2 × 10⁻⁴, because the leader side then forfeits the evaluations it needs for exploitation. The gap curve of HEF-BSO is flat-bottomed between forty and sixty percent, whereas the rival's improves monotonically until its own leader side starves, a signature of warm-start coordination amortizing follower cost. On the constraint-coupling axis, the leader gap of HEF-BSO stays within a factor of two for scales in [0.5, 2], degrading to 1.7 × 10⁻⁵ at the 4× tightening, whereas BLEAQ improves to 7.8 × 10⁻⁵ under the 0.25× loosening, approaching the proposed method, but collapses to 3.2 × 10⁻² at the 4× tightening, a 38-fold degradation that traces to curvature loss in its quadratic response model.

*The 40/60 row is the reference configuration and reproduces the suite-wide medians of Table 4. Constraint-coupling scales s in [0.5, 2] leave the HEF-BSO median leader gap within a factor of two of its default, while BLEAQ degrades thirty-eight-fold at s = 4.

Behind the starvation headline stands the warm-start archive of Section 4.5 rather than brute follower accuracy: Each cached response places the next inner solve close to the true optimum, so a thirty-percent allowance still lands most solves within tolerance, and the constraint-coupling plateau spares practitioners a delicate tuning task. Classical bilevel experience warns that the coupling structure is exactly where hierarchical solvers fail quietly [55], [56], and the two-axis sweep quantifies that warning. A solver that certifies feasibility at a thirty-percent follower share needs far fewer inner iterations than one requiring half the budget, which at the subject level of the EEG measurement chain buys shorter calibration sessions without a proportional loss of decoding accuracy. Where bilevel coordination pays, in summary, is the starved and tightened regime, precisely the regime that deployed hierarchical design problems inhabit.

Three families of perturbations act on the evaluation pipelines, following the lineage of uncertain-environment evolutionary optimization [57] and classical robustness studies in stochastic search [27]. The first family corrupts every function evaluation additively:

f(x,y) = f(x, y) + ε, ε ∼ (0, σ ), (13)

with σ expressed as a fraction of the empirical range of f on the instance and swept over the values 0.02, 0.05, 0.08, and 0.12; the lower-level objective f and, in a second experimental arm, the upper-level objective F are corrupted in separate experiments, with the worse arm defining the reported degradation. The second family, parameter perturbation, multiplies each of six control parameters, namely population size, cluster count, the endpoints of the step-size schedule, covariance decay, budget split, and shrinkage intensity, by independent factors uniform on [0.75, 1.25] around the reference configuration, representing ±25 percent misspecification by a practitioner. The third family, adversarial lower-level distortion, replaces the follower's true response with a corrupted one:

∗ ∗ ∗ ∗

yadv (x) = y (x) + δ (x), ‖δ ‖ ≤ Δ, (14)

where the adversary watches the leader distribution and chooses the distortion within a budget Δ = 0.15 in normalized lower-level coordinates so as to maximize the leader's optimality gap. The construction is an abstract model of systematic bias in the inner solver, such as an imperfect cross-validation split, a biased regularization path, or an artifact-contaminated fold. Throughout, degradation is reported as the relative increase of the median final gap with respect to the unperturbed control, so all three families share one reporting scale; the thirty runs per condition and instance preserve the statistical machinery of Section 6.4, and every algorithm receives the same perturbations identically.

Table 9 and Fig. 4 collect the outcomes. Under noise at σ = 0.08 of the function range, the median gap of HEF-BSO degrades by 9.4%, nested CMA-ES by 22.7%, BLEAQ by 26.3%, and KKT-SQP by 33.8%. The noise sweep of Fig. 4 shows a stable ordering from σ = 0.02 through σ = 0.12, with the degradation of HEF-BSO growing sub-linearly in σ, at 2.1 and 5.8% for the two smaller noise levels and 14.1% at the largest. Under ±25% parameter perturbation, the degradations are 7.6% for HEF-BSO against 17.8% for nested CMA-ES, 21.4% for BLEAQ, and 26.9% for KKT-SQP; for all population-based methods, the dominant sensitivity is the budget split, whose perturbation directly trades follower accuracy against leader exploration, followed by the covariance decay.

*Median final gap over 30 runs per condition on the ten SMD instances; degradation is the relative increase over the unperturbed control.

The degradation of HEF-BSO grows sub-linearly in σ and remains the lowest at every noise level.

Under adversarial lower-level distortion at Δ = 0.15, degradations reach 15.2% for HEF-BSO against 31.6% for nested CMA-ES, 38.9% for BLEAQ, and 47.2% for KKT-SQP. Doubly penalized is BLEAQ, because the distorted responses also enter the refits of its quadratic approximation, so the corruption compounds across generations, whereas a nested method merely inherits each distorted response as a biased evaluation. KKT-SQP remains the most fragile overall, since the corrupted response also feeds the complementarity residual that its SQP loop must absorb. Across all three families, the ratio of HEF-BSO's degradation to that of the strongest rival stays between roughly two-fifths and one-half, so the robustness advantage is of the same order as the accuracy advantage and is not an artifact of a single perturbation type.

The mechanism is best understood through variance scaling. Theory for evolution strategies under noisy evaluation establishes that the progress rate is maximized by a mutation strength that grows with the noise level, because larger steps average the noise over a broader region of the landscape [58]; covariance-aware mutation picks up this scaling automatically in the directions where the cluster's accepted ideas already disperse widely, and declines it where they do not. Shrinkage complements it: under noise the sample covariance of a small cluster is itself noisy, and the estimator of Eq. 8 trades variance for bias precisely in the regime where unbiased estimation is hopeless [36].

Eigenvalue clipping caps the damage that an unlucky noise-driven update can inflict on the mutation geometry, confining the condition number of every search distribution and thereby guaranteeing the non-degeneracy Assumption 3 of Section 5.1 by construction rather than by hope. Warm-start coordination adds a subtler buffer: A follower started near the true response turns a perturbed evaluation into a small correction rather than a full re-solve, which truncates the propagation of lower-level noise into the leader's incumbent. The net effect is that perturbations are absorbed as rescaled, not amplified, search variance, and this operational meaning of robustness, namely bounded degradation under an explicitly enumerated perturbation protocol, is proposed as a reporting standard for bilevel metaheuristics.