A Structural and Methodological Audit of Quantum-Inspired Differential Evolution

1. Introduction

In this paper, we analyze multiple quantum-inspired versions of Differential Evolution (DE) [1,2,3] and demonstrate that their published formulations contain structural ambiguities that prevent unambiguous implementation. In the absence of publicly available source code, independent replication relies entirely on the provided pseudocodes and equations. In Section 2, we illustrate that this process is hindered by mathematical inconsistencies and undefined parameters. Because these publications [1,2,3] share a foundational framework, these structural inconsistencies propagate across the proposed algorithms. We provide a visual summary of these inconsistencies in Figure 1, highlighting several instances where the mathematical formulation is insufficient to establish convergence.

The continuous optimization literature features numerous “quantum-inspired” metaheuristics, where concepts such as superposition or amplitude collapse serve as structural metaphors for probabilistic sampling and population diversity. Employing these metaphors within classical algorithms is an established practice. However, the frameworks under investigation explicitly cross the boundary from conceptual inspiration to claims of physical quantum advantage. The authors repeatedly assert that the algorithm can “balance the side effects of quantum computation” and “fully leverage the advantages” of quantum processing to mitigate high-dimensional complexity [1,2,3].

The first of these variants appears in the paper An Enhanced MSIQDE Algorithm With Novel Multiple Strategies for Global Optimization Problems,deng2020enhanced. There are several claims that “the quantum-inspired DE (QDE) fully utilizes the rapidity of quantum computing and the optimization performance of DE” and that “it can avoid premature, accelerate convergence, and guarantee the diversity of the population”.

A subsequent quantum-inspired framework is introduced in the paper Quantum differential evolution with cooperative coevolution framework and hybrid mutation strategy for large scale optimization,deng2021quantum. The authors utilize the cooperative coevolution (CC) framework to develop CCQDE. They also design a hybrid mutation strategy based on local neighborhood mutation and SaNSDE.

In a recent publication, Deng et al. proposed the PSEQADE [3] algorithm to solve large-scale global optimization problems. The authors assert that by "taking advantage of the characteristics of quantum computation," the algorithm mitigates dimensional scaling issues, explicitly claiming that its convergence performance is “hardly affected by a significant increase in dimensions” and effectively overcomes the “curse of dimensionality”.

We provide a detailed analysis of the most recent PSEQADE [3] algorithm. Our mathematical and algorithmic audit reveals that the claims about algorithms’ exceptional performance are challenged by a series of structural and methodological inconsistencies: contradictory equations, a divergent scale factor, undefined operators, deviations from standard benchmark protocols, and unusual statistical outcomes. The paper [3] proposes that PSEQADE addresses the poor optimization ability of traditional DE and its variants used in CEC competitions. That is a very strong statement, particularly when heavily adapted classical DE variants continue to demonstrate robust performance on complex IEEE CEC benchmark landscapes.

2. Mathematical and Algorithmic Flaws

The PSEQADE [3], EMSIQDE [1], and HMCFQDE [2] algorithms rest on a theoretical foundation undermined by structural inconsistencies. The following subsections detail multiple analytically provable mathematical failures within the published equations. It is critical to note that these are not minor inefficiencies or suboptimal parameter choices; rather, each inconsistency substantially complicates implementation and raises questions regarding reproducibility.

2.1. Structural Ambiguities, Dimensional Mismatches, and Variable Overloading

To evaluate the mathematical consistency of the proposed framework, it is necessary to examine the foundational variable upon which these quantum-inspired variants rely. The individuals are encoded as “quantum chromosomes” composed of “qubits,” mapping a latent representation into the real-world problem solution space $\mathsf{\Omega }=[a,b]$ [2,3].

This core transformation step, presented as Eq. (12) in [2,3] and Eq. (9) in [1], is formulated as follows:

$$q_{i,j}={\displaystyle \frac{(b-a)·q_{i,j}+(b+a)}{2}}={\displaystyle \frac{1}{2}}\left[\begin{array}{c}(b-a)·cos\left({\theta }_{i,j}\right)+(b+a)\\ (b-a)·sin\left({\theta }_{i,j}\right)+(b+a)\end{array}\right]$$

(1)

An analysis of this expression reveals several foundational inconsistencies from its initial definition, which may stem from compounding typographical omissions or structural oversights.

First, the expression introduces an unindexed dependency where the target variable $q_{i,j}$ appears simultaneously on both sides of the equality sign without an iterative or temporal subscript (such as t or $t+1$). If this formulation is treated as a static algebraic identity, isolating and solving for $q_{i,j}$ yields:

$$2q_{i,j}=(b-a)q_{i,j}+b+a\Rightarrow q_{i,j}={\displaystyle \frac{b+a}{2-b+a}}$$

(2)

This algebraic constraint reduces $q_{i,j}$ to a fixed scalar constant determined entirely by the domain boundaries, structurally decoupling it from the stochastic evolutionary process and the latent phase angle ${\theta }_{i,j}$.

Second, the equation introduces an immediate dimensional mismatch. The leftmost expression defines $q_{i,j}$ as a scalar or standard coordinate component, whereas the rightmost expression defines it as a $2\times 1$ column vector containing both cosine and sine branches. This conflation makes it mathematically ambiguous whether $q_{i,j}$ is intended to represent a spatial coordinate, a probability amplitude vector, or a dual-component state.

Beyond the initial formulation, the exact physical and mathematical definition of $q_{i,j}$ shifts arbitrarily across sequential operations, a practice that introduces significant symbol overloading. Across Equations (12), (13), and (22) in the manuscript of Deng et al. [3], the identifier $q_{i,j}$ is redefined to represent three distinct, mutually incompatible mathematical objects:

• Phase 1 (State Vector / Coordinate Mapping): In Equation (12), it acts as a $2\times 1$ multi-component object containing trigonometric mappings, while simultaneously referencing its own pre-transformation state on the right-hand side [3].
• Phase 2 (Scalar Probability Amplitude): In the subsequent observation step described by Equation (13), $q_{i,j}$ is processed strictly as a scalar probability amplitude [3].
• Phase 3 (Latent Phase Angle): In Equation (22), the same variable is utilized to denote a deterministic latent phase angle derived from an inverse operation [3].

This systematic conflation of a spatial coordinate, a multi-dimensional amplitude vector, and a phase angle within a single identifier removes dimensional consistency from the algorithmic description. Because the exact nature of this core variable changes depending on the execution step, the underlying transformation sequence remains unresolvable from the text alone.

2.2. Numerical Divergence of the Adaptive Scale Factor

In Section 3.1 of [3], the “quantum-adaptive mutation strategy” is introduced, intended to dynamically reduce mutation intensity. The scale factor F is defined in Equation (14) as

$$F=l{f}_0·e^{1+{\displaystyle \frac{G_m}{G_m-G+1}}},$$

where $G_m$ is the maximum number of iterations, G is the current iteration, and $l{f}_0=0.5$. There is an explicit claim that F starts at 0.5 and “ends at infinitely close to 0” to balance the effects of excessive mutation.

A monotonicity check of the exponent demonstrates that as the iteration count G increases, the denominator $(G_m-G+1)$ decreases toward unity

$${\displaystyle \frac{dF}{dG}}=l{f}_0·e^{1+{\displaystyle \frac{G_m}{G_m-G+1}}}·{\displaystyle \frac{G_m}{{(G_m-G+1)}^2}}>0.$$

Because the derivative is strictly positive, F is a strictly increasing function. This directly contradicts the claim that F “decreases with the increasing number of iterations”.

The implication of this is the numerical stability of the PSEQADE algorithm. Evaluating the limit as $G\to G_m$ yields

$$\underset{G\to G_m}{lim}F=l{f}_0·e^{1+G_m}=0.5·e^{1+G_m}.$$

For a standard run of $G_m=1000$, F starts in the first generation at $\approx 3.695$ and reaches a value of $\approx 0.5·e^{1001}\approx {10}^{434}$, which far exceeds the maximum representable value for a 64-bit double-precision float ($\approx 1.8\times {10}^{308}$). Far from converging to zero, the proposed scale factor triggers immediate floating-point overflow, which would theoretically render the operator non-functional without undocumented interventions.

2.3. Probabilistic Inconsistency of the Sampling Operator

The observation mechanism in [3] contains an internal logical inversion that is inconsistent with the paper’s own probabilistic definitions. In Section 2.2 of [3], the amplitude is defined as $\alpha =cos\left({\theta }_{ij}\right)$, so by standard Born-rule sampling, the cosine branch should be selected with probability ${cos}^2\left({\theta }_{ij}\right)$. Equation (13) of the original manuscript instead proposes:

$$q_{ij}=\left\{\begin{array}{cc}cos\left({\theta }_{ij}\right),\hfill & \mathrm{if}{cos}^2\left({\theta }_{ij}\right)<\mathrm{rand}(0,1)\hfill \\ sin\left({\theta }_{ij}\right),\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

Because $\mathrm{rand}(0,1)\sim U(0,1)$, the condition ${cos}^2\left({\theta }_{ij}\right)\mathrm{rand}(0,1)$ holds with probability $1-{cos}^2\left({\theta }_{ij}\right)={sin}^2\left({\theta }_{ij}\right)$. The operator therefore selects $cos\left({\theta }_{ij}\right)$ with probability ${sin}^2\left({\theta }_{ij}\right)$ and $sin\left({\theta }_{ij}\right)$ with probability ${cos}^2\left({\theta }_{ij}\right)$, the opposite of what the paper’s own amplitude definition requires. This self-contradiction is independent of any external sampling convention and renders the observation step internally inconsistent.

2.4. Initialization Range and Structural Positive Bias

The algorithm restricts the initial quantum phase angle to ${\theta }_{i,j}=\pi ·\mathrm{rand}(0,1)$ in Equation (11) [3], which limits the latent initialization space to $[0,\pi ]$. Because $sin\left(\theta \right)\ge 0$ across this entire interval, any step that selects the sine branch during observation yields an exclusively non-negative amplitude. When these amplitudes are subsequently mapped to a symmetric real-world solution space via the transformation equations, the population is structurally blocked from uniformly sampling the lower half of the domain at initialization.

The mathematical proof of this systemic asymmetry is revealed by evaluating the analytical expectation of the initialized amplitude q under a uniform distribution $\theta \sim \mathcal{U}(0,\pi )$. Incorporating the written observation probabilities, the expected value is derived as:

$$\mathbb{E}\left[q\right]={\displaystyle \frac{1}{\pi }}{\int }_0^{\pi }\left[{cos}^3\left(\theta \right)+{sin}^3\left(\theta \right)\right]d\theta ={\displaystyle \frac{4}{3\pi }}\approx 0.424$$

This non-zero expectation deviates from the ideal value of zero required for an unbiased, symmetric initialization. It structurally shifts approximately 75% of the initial sampled values into the positive quadrant of the search space, whereas a standard interval choice such as $[-\pi ,\pi ]$ or $[0,2\pi ]$ would naturally produce an unbiased, near-equal positive and negative split.

The immediate consequence of this restricted boundary choice is an artificial spatial bias that systematically biases the initial population toward the positive quadrant of the objective landscape. For an optimizer to successfully navigate symmetric or origin-centered benchmark functions from such an asymmetrical starting configuration, the downstream mutation operators must forcibly pull the population back toward zero.

2.5. Mathematical Realization and the Classical 2D Coordinate Representation

Evaluating the operational population update mechanism (distributed across Equations 11 and 12) reveals that these claims of computational advantage are not supported by the underlying mechanics. Rather than executing true quantum operators on a physical statevector, the formulation relies entirely on classical trigonometric parameters. In standard quantum information theory, a single-qubit state space operates within a two-dimensional complex Hilbert space ${\mathbb{C}}^2$. The associated complex amplitudes encapsulate both a probability magnitude and a relative phase factor, providing two degrees of freedom that map to a three-dimensional Bloch sphere. Consequently, state transformations are governed by unitary operators from the special unitary group $SU\left(2\right)$.

The framework under review restricts its state parameters strictly to real numbers ($\mathbb{R}$). By defining the state vector as the column vector ${[cos\theta ,sin\theta ]}^T$, the relative phase vector and its associated degrees of freedom are eliminated. Because the algebraic intersection of the special unitary group with real matrices yields the special orthogonal group ($SU\left(2\right)\cap {\mathbb{R}}^{2\times 2}=SO\left(2\right)$), the core operator structurally reduces from a quantum Bloch sphere transformation to a classical two-dimensional coordinate rotation matrix. This mathematical reduction demonstrates that the algorithm does not employ quantum computational primitives, rendering assertions of intrinsic quantum advantage or quantum-specific scaling benefits structurally unsupported.

2.6. Operational Formulations of the Rotation Mechanism and Missing Parameters

The practical implementation of the core search operator reveals a substantial contradiction between the high-level claims of balancing global and local search and the actual mathematical details provided in the text. In [1], the visual and mathematical presentation of the state update sequence contains significant structural ambiguities. The mandatory Dirac notation (〉) is systematically replaced with a greater-than-or-equal-to inequality operator (≥), rendering the state vectors as $|{\phi }^{\prime }\ge$ and $|\phi \ge$. Furthermore, evaluating the linear algebraic transformation reveals that the matrix product layout in the text is structurally disjointed; Equation (12) in [1] does not constitute a self-contained equation, but rather an unaligned fragment of the trigonometric expansion from the preceding row.

Figure 2. Reproduction of the state update layout from [1]. The notation utilizes inequality symbols in place of standard right-hand Dirac brackets, and the secondary line contains an isolated matrix fragment resulting from layout overflow. The underlying operation reduces a purported quantum primitive to a classical two-dimensional coordinate rotation.

Figure 2. Reproduction of the state update layout from [1]. The notation utilizes inequality symbols in place of standard right-hand Dirac brackets, and the secondary line contains an isolated matrix fragment resulting from layout overflow. The underlying operation reduces a purported quantum primitive to a classical two-dimensional coordinate rotation.

Expanding the complete matrix product reveals that the operation reduces entirely to a classical, scalar addition of the phase angle ($\theta \leftarrow \theta +\Delta \theta$). Characterizing this straightforward trigonometric shift as an advanced quantum mechanism introduces a misleading presentation of the algorithm’s actual computational primitives, masking a purely classical coordinate adjustment behind quantum nomenclature.

The rotation mechanism is cited across multiple publications [1,2,3] as a primary driver of convergence and exploratory balance. For instance, it is explicitly invoked in Step 10 of the main loop in [3] to balance global and local search dynamics. However, across all three manuscripts, the parameter $\Delta \theta$ is introduced solely through a duplicated conceptual schematic figure depicting a phase rotation. No mathematical formula, adaptive update rule, or concrete initialization value for calculating or bounding the magnitude of $\Delta \theta$ is provided anywhere in the text, pseudocode, or figure captions.

3. Empirical Discrepancies and Methodological Limitations

Implausible Convergence and Zero Variance

In Tables 6 through 10 of the original paper [3], the authors report that PSEQADE achieves a "Best", "Mean", and "Std" of exactly $0.00E+00$ on numerous complex landscapes (e.g., $f_1,f_4,f_5,f_8,f_9$) at 500, 1000, and 3000 dimensions.

Shifting: 

In the actual CEC 2017 suite [4], shifting ensures that the global optimum is not located at the origin $(0,0,...,0)$, preventing algorithms with inherent “center-seeking” biases from artificially inflating performance. In contrast, Table 1 of the manuscript lists $f_{min}=0$ for all 20 evaluated functions. This further confirms that the authors did not utilize the official CEC 2017 benchmark code or specifications, but rather evaluated their algorithm on basic functions centered at zero.

Omission of Coordinate Rotation and Non-Separability 

The experimental validation is significantly compromised by the removal of coordinate rotation from the CEC 2017 benchmark suite. In continuous optimization, rotation matrices are strictly applied to base functions to induce non-separability among variables. This non-separability is a defining hurdle of modern test landscapes, explicitly designed to prevent naive algorithms from succeeding via simple dimension-by-dimension coordinate descent. To effectively navigate these correlated parameter interactions, genuine state-of-the-art optimizers are forced to utilize sophisticated mechanisms, such as eigen crossover or covariance matrix adaptation.

Function evaluation starvation of baselines. 

The comparison algorithms, MLSHADE-SPA, SHADE-ILS, and CCPSO2, are large-scale global optimisation algorithms designed for evaluation budgets on the order of $3\times {10}^6$ function evaluations on 1000-dimensional problems. The authors provide 100,000 function evaluations for 3000-dimensional problems, equivalent to approximately 33 evaluations per dimension. At this budget, these algorithms have not completed a single meaningful adaptation cycle.

In Table 1 (here in appendix Table A1) of the evaluated manuscript, the authors present the standard, unshifted, and unrotated mathematical formulas for the base functions. For example, their definition of $f_1\left(x\right)=x_1^2+{10}^6{\sum }_{i=2}^D{x}_i^2$ contains no shift vector or rotation matrix, leaving the global optimum explicitly at the origin.

Further Simplification of Hybrid Functions: 

In the official CEC 2017 suite, Hybrid Functions (e.g., $F_{10}$ through $F_{19}$) are constructed by randomly dividing the variables into subcomponents and applying different basic functions to each subcomponent. The authors define their hybrid functions ($f_{15}$ through $f_{20}$) as simple, global weighted sums of basic functions (e.g., $f_{15}\left(x\right)=0.2f_2\left(x\right)+0.4f_3\left(x\right)+0.4f_4\left(x\right)$). This formulation fails to test the algorithm’s ability to handle subcomponents with varying properties.

3.1. Minor Algorithmic and Parametric Ambiguities

Beyond the primary structural inconsistencies, the text contains several secondary omissions that further obscure its implementation. Regarding parameterization, the crossover rate ($CR$) is entirely omitted from the parameter tables and text despite its known sensitivity in high-dimensional DE settings, and critical Population State Evaluation thresholds ($D{T}_{\mathrm{init}}$ and $A{T}_{\mathrm{init}}$) are either completely unstated or uncalibrated for the specified domains. Structurally, the complexity analysis deviates from standard asymptotic notation by retaining lower-order constants, such as $O\left((O_f\left(D\right)+6D+3)NP\right)$, while simultaneously omitting the mandatory generation multiplier ($G_{max}$) required to express total runtime complexity rather than a single-generation cost. Finally, the empirical parameter configurations feature an unaddressed asymmetry, halving the total iteration allowance for baseline competitors relative to PSEQADE. While these details do not alter the core mathematical contradictions of the operators, their omission or misconfiguration reinforces the unresolvable nature of the published framework.

4. Reproducibility Issues

The independent reproduction was implemented in Python using a population of $NP=50$ individuals and a crossover rate of $CR=0.9$, with the maximum number of generations derived as $G_{max}=\lfloor \mathrm{MaxFEs}/NP\rfloor$ to match the original evaluation budget exactly. The PSE execution period was set to $IT=250$ generations, with the aggregation threshold initialised as $A{T}_{\mathrm{init}}=0.02$ and the diversity threshold derived as $D{T}_{\mathrm{init}}=IT/10000=0.025$; none of these values are stated in [3] and represent necessary assumptions.

Due to the absence of source code and the structural ambiguities documented in the preceding sections, several implementation decisions required explicit correction or assumption:

1.

Scale factor: The published Eq. (14) diverges to ${10}^{434}$ under standard execution, the implementation substitutes $F=l_{f0}·exp(-G/(G_{max}-G+1))$, the only formulation consistent with the paper’s stated qualitative intent of smooth decay toward zero, and presumed to reflect a sign typographic error in the exponent.

2.

Quantum rotation gate: No formula for the step-size $\Delta \theta$ is provided anywhere in [3]; the implementation defines $\Delta {\theta }_{i,j}=s·GP·({\theta }_{\mathrm{best},j}-{\theta }_{i,j})$ where $s=0.05$ is an assumed rotation strength and $GP=(G_{max}-G+1)/G_{max}$ is a linearly decaying progress ratio, chosen to reflect the paper’s description of balancing global and local search over time.

3.

Observation operator: The inverted inequality in Eq. (13) was corrected to $\mathrm{rand}(0,1)<{cos}^2\left({\theta }_{ij}\right)$ to restore consistency with the paper’s own amplitude definitions.

4.

PSE thresholds: $D{T}_{\mathrm{init}}$ and $A{T}_{\mathrm{init}}$ are uncalibrated in the original manuscript and were set heuristically.

Random seeds were fixed via numpy.random.seed before each independent run. The full implementation is available at https://github.com/VojtechNovak/PSEQADE.

We showcase convergence in Figure 3 using multiple CEC-winning optimizers, SciPy implementation of DE, and iL-SHADE from the pyade package. Every optimizer converges in high precision, but PSEQADE does not.

5. Conclusions

This paper has presented a systematic mathematical and methodological audit of the PSEQADE algorithm proposed by Deng et al. [3]. The analysis has outlined several inconsistencies in the mathematical formulation, distinct divergences in the benchmarking methodology, and a persistent contradiction between the high-level claims of quantum utility and the purely classical, trigonometric nature of the underlying operators. Taken together, these findings indicate that the published formulation is insufficient to establish the algorithm’s stated functionality or to permit unambiguous replication.

The mathematical foundation of the framework is characterized by structural ambiguities across its primary operators. First, the initialization range introduces an inherent positive spatial bias that restricts uniform sampling. Second, the observation mechanism features an internal logical inversion, rendering its behavior inconsistent with standard amplitude-based sampling or its own probabilistic interpretation. Third, the foundational variable $q_{i,j}$ is subjected to significant symbol overloading, shifting between a spatial coordinate, a multi-dimensional vector, and a latent phase angle in a manner that removes dimensional consistency. Furthermore, the adaptive scale factor F possesses a strictly positive derivative, causing it to monotonically increase and diverge toward ${10}^{434}$, which triggers immediate floating-point overflow under standard double-precision execution. Finally, the central rotation gate lacks an operational formula or step-size update rule, while the population state evaluation thresholds remain structurally uncalibrated for the intended experimental context.

During empirical validation, baseline competitors were subjected to evaluation starvation. Additionally, the mathematical structure of the metrics used strongly implies that the evaluation was conducted on simplified, unshifted, and unrotated base functions rather than the official, highly correlated benchmark configurations.

The global optimization community has spent decades refining highly sophisticated metaheuristics, with state-of-the-art adaptive Differential Evolution variants like L-SRTDE [5] established through rigorous IEEE CEC competitions representing the current performance frontier. Results that comprehensively surpass these baselines on complex, multi-modal, and high-dimensional landscapes would therefore represent an extraordinary advancement requiring correspondingly rigorous methodological support. The reported results, an absolute mean fitness of $0.00\mathrm{E}+00$ and standard deviation of $0.00\mathrm{E}+00$ across 25 independent runs under a budget of only $50,000$ function evaluations [3], cannot be reconciled with the published formulation as stated.