1
A Family of Quaternion-Valued Differential Evolution Algorithms for Numerical Function Optimization
arXiv:2605.12362v1 [cs.NE] 12 May 2026
Gerardo Altamirano-Gomez1 , Álvaro Gallardo2 and Carlos Ignacio Hernández Castellanos1
Abstract—The numerical optimization of continuous functions is a fundamental task in many scientific and engineering domains, ranging from mechanical design to training of artificial intelligence models. Among the most effective and widely used algorithms for this purpose is Differential Evolution (DE), known for its simplicity and strong performance. Recent research has shown that adapting AI models to operate over alternative number systems-such as complex numbers, quaternions, and geometric algebras-can improve model compactness and accuracy. However, such extensions remain underexplored in bio-inspired optimization algorithms. In particular, the use of quaternion algebra represents an emerging area in computational intelligence. This paper introduces a family of novel QuaternionValued Differential Evolution (QDE) algorithms that operate directly in the quaternion space. We propose several mutation strategies specifically designed to exploit the algebraic and geometric properties of quaternions. Results show that our QDE variants achieve faster convergence and superior performance on several function classes in the BBOB benchmark compared to the traditional real-valued DE algorithm. Index Terms—Quaternions, Differential Evolution, Evolutionary Computation, Numerical function optimization, Representation of individuals
T
HE last decades have seen an explosion in the creation of bio-inspired algorithms for solving optimization problems. The task of designing such algorithms involves finding a suitable representation for the problem and developing an appropriate search method. For numerical optimization of continuous N -dimensional functions, a N -dimensional vector is the most common and natural way of representing the solution; however, in the last years, the use of alternative representations based on complex numbers [1] or quaternions [2]–[6] has been explored. These works can be classified in two categories, according to how they apply the quaternion representation: In the first category, the quaternion representation is used for modeling the problem, but the solution of the optimization problem is obtained with the real-valued versions of the algorithms. For example, Song et al. [6] use this methodology for modeling the protein-ligand docking problem, while Saldias et al. [3] model a human knee joint. Both works obtained the optimal solution using the real-valued DE algorithm. • In the second category, the algorithm associates each value of an N -dimensional vector with a quaternion; •
1 Instituto de Investigaciones en Matemáticas Aplicadas y Sistemas, Universidad Nacional Autónoma de México, México 2 Universidad Iberoamericana, México
it finds the minimum in a 4-dimensional space using the real-valued versions of the optimization algorithms, and obtain the solution to the problem by mapping each quaternion to a real number. This methodology has been applied in the Quaternion-valued versions of the Genetic algorithm [5], the Firefly algorithm [2], and the Bat algorithm [4]. Despite the advances in this direction, the task of designing a bio-inspired algorithm with an appropriate search method that works in the quaternion space remains an open problem. Moreover, other artificial intelligence algorithms, such as Support Vector Machines [7], Multi-Layer Perceptrons [8], Convolutional Neural Networks [9], [10], and other Neural Networks models [11]–[13] have exploited the properties of the quaternion representation by adapting their inner machinery to the use of quaternion algebra. In particular, it has been shown that Quaternion-valued CNNs converge faster during training due to the “navigation on a much compact parameter space during learning” [14], and for Quaternion SVMs, it was proved that splitting the hypercomplex numbers into parts, and then solving the primal and dual problems independently in the real domain, loses the benefits that the hypercomplex numbers have when are embedded in their spaces, e.g. sparsity, and taking advantage of geometric and topological features [7]. Importantly, the use of quaternion representations in optimization and learning algorithms constitutes an emerging area within computational intelligence. This growing interest is driven by their algebraic richness, compactness, and ability to represent multi-dimensional transformations naturally. To address this gap, we propose a Quaternion-Valued Differential Evolution Algorithm (QDE), which search method is embedded in the quaternion space and exploits the properties of the quaternion algebra. The main contributions of this work are as follows: We propose a family of QDE algorithms that operate entirely in quaternion space. • We propose two initialization mechanisms: Euclidean and polar. • We propose six novel mutation strategies that leverage the algebraic and geometric properties of quaternions, including both Euclidean and polar-based mechanisms. • We conduct a statistical analysis using the BBOB benchmarksuite, identify the most effective mutation strategies and initialization methods, and demonstrate that our QDE algorithms achieve faster convergence and comparable
•
2
Fig. 1. Two similar triangles, turned and in a common plane. a) General configuration b) Right triangles with equal length cathetus configuration. Adapted from [16].
or superior performance relative to the classical DE algorithm. The rest of the paper is organized as follows. The foundations of Quaternion Algebra is presented in Section I followed by the Quaternion-Valued Differential Evolution Algorithm in Section II. Then, the Experimental Analysis and the Discussion of Results are presented in Sections III and IV, respectively. Finally, Section V states the Conclusions and Future Works. I. Q UATERNION A LGEBRA This mathematical system was developed by W.R. Hamilton (1805-1865) at the middle of the XIX century [15]–[17]. This section presents A brief introduction to quaternions as was developed by W.R. Hamilton. His work on this subject started by exploring ratios between geometric elements, consequently he called quaternion to the quotient of two vectors. Next, we summarize the construction of quaternions as was introduced by Hamilton [16]. Let AOB and COD be two similar triangles lying on a common plane, which are similarly turned, ~ OB, ~ OC, ~ and OD ~ be the vectors from point O and let OA, to point A, B, C, and D, respectively, see Figure 1a); then, the ratio of the vectors satisfy the following equality [16, pp. 112]: ~ : OA ~ = OD ~ : OC, ~ OB (1) which can be expressed as a geometric fraction: ~ ~ OD OB = . ~ ~ OA OC
(2)
Then, by making the triangle COD into BOA′ , see Figure 1b), the following relationship holds [16, pp. 130]: ~ ′ ~ OA OB = . ~ ~ OA OB
(3)
~ OA, ~ Thereafter, we multiply each side of the equation by OB/ thus: !2 ~ ′ ~ OA OB = . (4) ~ ~ OA OA ~ and OA ~ ′ have the same magnitude, but opposite Since OA ~ = −OA ~ ′ , then: direction, that is OA !2 ~ OB = −1. (5) ~ OA
~ and OB, ~ we construct the following geometric Fig. 2. Given two vectors OA ~ ′′ lies on the plane ~ vector OB configuration: plane Π is orthogonal to OA, ~ ′ is the projection ~ into Π, while vector OB OAB and is the projection of OB ~ into OA. ~ Adapted from [16]. of OB
Equation 5 shows that the quotient of two perpendicular vectors of equal length, that is a quaternion, equals the square roots of negative unity [16, pp. 131]. In addition, for a quaternion: q = v~1 : v~2 , where v~1 and v~2 are vectors, we can rewrite it as qv~2 = v~1 . In this case, q is called a versor, i.e. an element that transforms v~2 into v~1 by rotating it [16, pp. 133]. In this way, Hamilton connected the quaternion representation √ with the root −1, and its geometric meaning. Now, we are going to connect the quaternion concept with ~ : its modern representation. Considering Figure 2, let q = OB ′ ′′ ~ ~ ~ OA be a quaternion, and let OB and OB be parallel and ~ respectively; such that OB ~ = perpendicular vectors to OA, ′ ′′ ~ ~ OB + OB . Then, we can decompose the quaternion, q, into: ~ ′ : OA ~ + OB ~ ′′ : OA. ~ q = OB
(6)
~ ′ and OA ~ are parallel vectors, their quotient is just Since OB ~ into OB ~ ′, a scale factor representing the projection of OA so first term turns into a scalar. The second term represents ~ on the plane through O, which is the projection of OA ~ in addition, this means OB ~ ′′ can be perpendicular to OA; ~ by applying a versor transformation. Since obtained from OA ~ ′′ and OA ~ are perpendicular to each other, they can be OB expressed as a linear combination of right versors (unitary quaternions, orthogonal to each other). This leads to the wellknown result that every quaternion can be expressed in the quadrinomial form [16, pp. 160]: q = qR + qI î + qJ ĵ + qK k̂,
(7)
where qR , qI , qJ , qK are scalars, and î, ĵ, k̂ are orthogonal unit quaternions. In terms of modern mathematics, the quaternion algebra, H, is: the 4-dimensional vector space over the field of the real numbers, generated by the basis {1, î, ĵ, k̂}, and endowed with
3
the following multiplication rules (Hamilton product): (1)(1) = 1,
TABLE I D EFINITION OF MATHEMATICAL SYMBOLS .
(8) Symbol R H bold lowercase letter, e,g. q qR
(1)(î) =ĵ k̂ =−k̂ ĵ = î, (1)(ĵ) = k̂ î = −îk̂ = ĵ,
(1)(k̂) = îĵ = −ĵ î = k̂,
î2 = ĵ 2 = k̂ 2 =−1.
(9)
The quaternion algebra is associative and non-commutative, and for two arbitrary quaternions: p = pR + pI î + pJ ĵ + pK k̂ and q = qR +qI î+qJ ĵ+qK k̂, their multiplication is calculated as follows: pq =pR qR − pI qI − pJ qJ − pK qK +
(pR qI + pI qR + pJ qK − pK qJ )î+
(pR qJ − pI qK + pJ qR + pK qI )ĵ+ (pR qK + pI qJ − pJ qI + pK qR )k̂.
(10)
Notice that each coefficient of the resulting quaternion, is composed of real and imaginary parts of the factors p and q. In this way, the Hamilton product capture inter-component relationships between both factors. Next, there are introduced some useful operations with quaternions. Let, q = qR +qI î+qJ ĵ+qK k̂, be a quaternion, its conjugate is defined as: q̄ = qR − qI î − qJ ĵ − qK k̂.
(11)
And its magnitude is computed as follows: p kqk = qq̄.
(12)
As well as the complex numbers, quaternions can be represented in polar form [18], [19], as follows: " # qI î + qJ ĵ + qK k̂ qθ = kqk cos(θ) + sin(θ) , (13) kqI î + qJ ĵ + qK k̂k
where: θ = atan
p
2 qI2 + qJ2 + qK qR
!
.
(15)
where wθ is a quaternion expressed in polar form: wθ = cos(θ) + sin(θ)(wI î + wJ ĵ + wK k̂),
(16)
and applies a rotation, with angle θ, along the axis wI î + wJ ĵ + wK k̂. Alternatively, we can split the transformation as a sandwiching product: p = w θ qw̄ θ ; 2
2
î, ĵ, k̂ ˆ· ~· ¯·
which is isomorphic to the 4D rotation group SO(4) [19]. In addition, there exist a two to one homomorphism with the rotation group SO(3) [20], [21]. Finally, Table I summarizes the notation that will be used in the rest of this paper.
II. Q UATERNION - VALUED D IFFERENTIAL E VOLUTION A LGORITHM (QDE) This section presents a novel algorithm that applies the quaternion algebra for representing the solutions of an optimization problem, as well as for searching the solution in the quaternion space. The algorithm is presented for dealing with problems that have four components. This can be a 4-dimensional problem, but also, due to the isomorphism between R3 and the imaginary part of the quaternion, it can be applied to 3-dimensional functions if we set the real part of the quaternion to zero. In addition, the algorithm can be extended for problems with a dimension that is multiple of four, provided the solution is divided in groups of four dimensions, i.e. if ~x is the solution vector of a D-dimensional problem, it can be divided in D/4 quaternions [q1 , q2 , . . . , qD/4 ]; thereafter, we can simultaneously apply the following steps for each quaternion.
(14)
A quaternion can represent a geometric transformation, which is applied as follows: p = wθ q,
qI , qJ , qK
Meaning The filed of the real numbers The quaternion algebra A quaternion The real component of a quaternion q The imaginary components of a quaternion q The imaginary bases of H An unitary vector A vector Conjugation sign
(17)
in this case, the angle of each quaternion, wθ/2 , is divided to half. From the group theory perspective, the set of unitary versors lies on a 3-Sphere, S3 , embedded in a 4D Euclidean space [20]; and together with the Hamilton product form a group,
A. Initialization The first step consists of creating a set of quaternions. For this, we propose two types of initialization: 1) E4. A random number generator with uniform distribution assigns a real number to each element of the 4dimensional vector. 2) Polar. In this case, we apply the polar representation of the quaternion, see Equation 13. A random number generator assigns a real value between [−2π, 2π] for the angle, and a randomly generated unitary vector for the direction of the quaternion. In both cases, the generator follows a uniform distribution. The result is a unitary quaternion. Thereafter, we iterate between the process of mutation, recombination and selection for creating new populations until the optimum quaternion is located.
4
B. Mutation Once the population has been initialized, we apply a differential mutation function to produce a quaternion. In this work, different mutation functions are proposed and evaluated. These can be divided in three types: the first one is a straightforward extension of the traditional mutation formula of the DE algorithm to the quaternion domain; the second type is an analogy of the traditional formula, but using the polar representation of the quaternions; and the third one is a simple formula for computing mutations in quaternion space. 1) Euclidean Mutation: Let q1 , q2 and q3 be three randomly chosen quaternions taken from the current population; then, the mutant quaternion, v, is computed using one of the following methods: •
Sum of Differences (ESD). This equation is an extension of the traditional formula, in which vectors have been replaced by quaternions, as follows: v = q3 + α(q2 − q1 ), α ∈ R.
(18)
The algorithmic version is shown in Algorithm 1. • Generalized Sum of Difference (EGSD). This is a generalization of the traditional formula for the quaternion space; in this case, the scale is not a real number but a randomly generated quaternion. This modification transforms the magnitude and the orientation of the solution in the quaternion space, since the Hamilton product is applied when we multiply the random quaternion by the expression (q2 − q1 ), as follows: v = q3 + qr (q2 − q1 ), qr ∈ H.
(19)
Note that when the imaginary part of the random quaternion is the zero vector, we obtained the previously presented formula; in addition, when all the elements are real numbers, we obtained the traditional formula of the real-valued DE algorithm. The algorithmic version is shown in Algorithm 2. 2) Polar Mutation: In the traditional mutation formula of the DE algorithm, we randomly select two vectors from the population and compute its difference; then, the difference is scaled. If we have a polar representation of a quaternion, an analogous procedure is to randomly select two quaternions from the population, q1 and q2 , and to compute the quaternion that transforms one into the other; thereafter, we can scale the magnitude and phase of the quaternion as follows: qr = α[cos(βθ) + sin(βθ)n̂], (20) q where θ = 2 arccos 1+c 2 , c is the dot product between the imaginary components of q1 and q2 , n̂ is a unit vector resulting of the cross product between the imaginary parts of quaternions q1 and q2 , and α, β ∈ R are the scale factors. In addition, the traditional mutation formula adds the scaled difference to a vector. If the scaled difference is added to one of the vectors used in its computation, the mutation promotes exploitation. In contrast, if it is added to a different vector, the scheme favors exploration. In the quaternion space, the
exploitation scheme is obtained when we apply qr as a versor transformation on one of the selected quaternions, q1 : v = qr q1 q̄r .
(21)
This scheme is named PM1, and its algorithmic version is shown in Algorithm 3. Alternatively we can apply an exploration scheme, which involves randomly selecting a different quaternion from the population, q3 , and applying qr as a versor transformation: v = qr q3 q̄r .
(22)
This scheme is named PM3, and the algorithmic version is shown in Algorithm 4. A combination of exploration and exploitation is obtained when we apply qr as a versor transformation on q1 , and then we sum q3 to the result: v = q3 + qr q1 q̄r .
(23)
This scheme is named PM13, and the algorithmic version is shown in Algorithm 5. 3) Simple Mutation (RQ): Since any coordinate of the quaternion space can be reached just using rotations, we propose a simple method to explore the performance of mutation when we apply this property. In this method, we generate a random unit quaternion, qr , and apply the sandwiching product on the selected quaternion, q1 , as follows: v = qr q1 q̄r .
(24)
This causes a random rotation on q1 in the quaternion space. The algorithmic version is shown in Algorithm 6. C. Quaternion Uniform Crossover Suppose that we are working with a D-dimensional problem, let ~x be a target vector took from the population, it can be divided into D/4 quaternions [q1 , q2 , . . . , qM/4 ]. In addition, let ~v be a mutant vector, it can be divided into D/4 quaternions [v1 , v2 , . . . , vD/4 ]. Then, the crossover builds a new solution by selecting some quaternions from the solution vector and others from the mutant vector, as follows: ( vi if (rand(0, 1) ≤ Cr) ui = , (25) qi otherwise where Cr ∈ [0, 1] is the crossover probability and 1 < i < D/4. The result is called a trial vector. Note that for problems that use a single quaternion, the crossover procedure consists in selecting between the mutant quaternion, v, and the target quaternion, q, and no exchange of information is involved. Thus, in this case this step is unnecessary. D. Selection This step consists in computing the fitness value of the trial and target quaternions, and selecting that with a lower value. Finally, the full version of the Quaternion-Valued DE algorithm is shown in Algorithm 7.
5
Algorithm 1 Euclidean Sum of Difference Mutation Algorithm (ESD) Require: Three quaternions, q1 , q2 and q3 ; and a scalar, α. Ensure: The mutated quaternion ,v. 1: v = q3 + α(q2 − q1 ) Algorithm 2 Euclidean Generalized Sum of Difference Mutation Algorithm (EGSD) Require: Three quaternions, q1 , q2 and q3 . Ensure: The mutated quaternion, v. 1: qr ← Generate a random 4D vector with uniform distribution. 2: v = q3 + qr (q2 − q1 ) Algorithm 3 Polar Mutation 1 Algorithm (PM1) Require: Two quaternions, q1 and q2 ; and two scalars, α and β. Ensure: The mutated quaternion, v. 1: c ← (q1I î + qq 1J ĵ + q1K k̂) · (q2I î + q2J ĵ + q2K k̂) 1+c 2: θ ← 2 arccos 2 3: n̂ ← (q1I î + q1J ĵ + q1K k̂) × (q2I î + q2J ĵ + q2K k̂)
4: qr ← α[cos(βθ) + sin(βθ)n̂] 5: v ← qr q1 q̄r
Algorithm 4 Polar Mutation 3 Algorithm (PM3) Require: Three quaternions, q1 , q2 and q3 ; and two scalars, α and β. Ensure: The mutated quaternion, v. 1: c ← (q1I î + qq 1J ĵ + q1K k̂) · (q2I î + q2J ĵ + q2K k̂) 1+c 2: θ ← 2 arccos 2 3: n̂ ← (q1I î + q1J ĵ + q1K k̂) × (q2I î + q2J ĵ + q2K k̂)
4: qr ← α[cos(βθ) + sin(βθ)n̂] 5: v ← qr q3 q̄r
Algorithm 5 Polar Mutation 1-3 Algorithm (PM13) Require: Three quaternions, q1 , q2 and q3 ; and two scalars, α and β. Ensure: The mutated quaternion, v. 1: c ← (q1I î + qq 1J ĵ + q1K k̂) · (q2I î + q2J ĵ + q2K k̂) 1+c 2: θ ← 2 arccos 2 3: n̂ ← (q1I î + q1J ĵ + q1K k̂) × (q2I î + q2J ĵ + q2K k̂) 4: qr ← α[cos(βθ) + sin(βθ)n̂] 5: v ← q3 + qr q1 q̄r
Algorithm 6 Simple Mutation Algorithm (RQ) Require: A quaternion, q1 Ensure: The mutated quaternion, v. 1: qr ← Generate a random 4D vector with uniform distribution. 2: v ← qr q1 q̄r
Algorithm 7 Quaternion-valued Differential Evolution Require: Population size N p, dimension D, crossover rate Cr. Ensure: Optimal solution. Np 1: Initialize population {[q1,i , . . . , qD/4,i ]}i=1 using one of: 2: E4: Random 4D vector generated with uniform distribution. 3: Polar: Generate a random angle θ ∈ [−2π, 2π] and a random unit vector direction with uniform distribution; then compute Equation 13. 4: repeat 5: for each set of quaternions [q1,i , . . . , qD/4,i ] in the population do 6: Generate distinct random indices: r0, r1, r2, such that: r0 6= r1 6= r2 6= i 7: jrand ← rand(1, D/4) 8: for j = 1 to D/4 do 9: if rand(0, 1) ≤ Cr or j = jrand then 10: uj,i ← M utation(qj,r0 , qj,r1 , qj,r2 ) 11: else 12: uj,i ← qj,i 13: 14: 15:
for i = 1 to N p do if f (u1,i , . . . , uD/4,i ) ≤ f (q1,i , . . . , qD/4,i ) then [q1,i , . . . , qD/4,i ] ← [u1,i , . . . , uD/4,i ]
16: until termination criterion is met
III. E XPERIMENTAL ANALYSIS This section presents the evaluation of the QDE algorithm and its comparison versus the real-valued DE algorithm. Thus, we empirically compared 12 different QDE schemes (the combinations of the two initialization methods with the six mutation functions) over a test-suit of 24 problems from the BBOB Benchmark [22]. This divides the test functions in five groups, where only the first one contains separable functions: 1) Separable functions. They can be decomposed into a sum or product of functions. Examples included in the benchmark are: the Sphere, the Ellipsoidal, the Rastrigin, the Büche-Rastrigin and the Linear Slope functions. 2) Functions with low or moderate conditioning and unimodal (U-Low). In this case, small changes in the input produce small changes in the output; thus, the function responds uniformly in all directions. Examples included in the benchmark are: the Attractive Sector, the Step Ellipsoidal, the Rosenbrock and the Rotated Rosenbrock functions. 3) Functions with high conditioning and unimodal (UHigh). These functions have a single optimal value, and the output of the function is sensitive to small perturbations in the input. Examples included in the benchmark are: the Ellipsoidal, the Discus, the Bent Cigar, the Sharp Ridge and the Different Powers functions. 4) Multi-modal functions with adequate global structure (M-Adequate). These functions have multiple local optima, but its landscape presents global clues, such as, regularities or symmetries; thus, the algorithms can navigate toward the global optimum. Examples included
6
TABLE II M EAN AND STANDARD DEVIATION OF THE FITNESS VALUES OBTAINED BY EACH ALGORITHM AT GENERATION 100 FOR S EPARABLE FUNCTIONS . Alg./Fn. E4-PM13 E4-ESD E4-PM1 E4-PM3 E4-RQ E4-EGS Polar-PM13 Polar-ESD Polar-PM1 Polar-PM3 Polar-RQ Polar-EGSD Real-DE
Rastrigin Median σ 10.311 4.358 7.481 3.022 1.298e-19 2.239e-21 1.297e-19 1.255e-21 7.637 4.853 10.983 4.809 2.638 2.798 2.515 5.372 1.322e-19 2.027e-21 1.324e-19 1.458e-21 2.013 1.022 2.019 4.927 0.995 0.671
Sphere Median σ 4.664e-19 1.952e-16 6.874e-10 1.727e-04 3.221e-39 2.304e-38 2.084e-39 1.511e-39 0.448 0.316 3.539e-03 5.583e-03 9.812e-25 1.555e-18 1.678e-09 1.44e-04 5.985e-40 1.934e-39 8.155e-41 3.602e-40 0.0346 0.04481 2.833e-04 4.751e-04 2.949e-11 8.223e-11
Ellipsoidal Median σ 0.384 20.723 9.378e-07 5.188e-03 8.559e-15 7.734e-12 8.563e-15 4.982e-18 29.544 85.428 226.527 180.791 0.01648 2.527 2.748e-08 0.01197 8.559e-15 9.417e-16 8.56e-15 4.992e-18 3.66 5.565 37.08 282.315 1.243e-07 4.665e-06
Buche-Rastrigin Median σ 2.901 3.016 3.115 1.862 9.539e-19 1.412e-20 9.633e-19 2.43e-20 9.194 3.551 10.36 5.204 8.875e-05 0.643 2.213 2.838 9.546e-19 2.076e-20 9.634e-19 1.041e-20 2.668 2.79 7.492 2.659 0.995 1.382
Linear Slope Median σ 0 0 0 5.882 57.686 4.892 53 6.698 48.79 6.183 0 0 1.46 8.053 0 5.081 61.614 1.198 61.749 1.771 60.563 0.556 0 0 0 0
TABLE III M EAN AND STANDARD DEVIATION OF THE FITNESS VALUES OBTAINED BY EACH ALGORITHM AT GENERATION 100 FOR FUNCTIONS WITH LOW OR MODERATE CONDITIONING . Alg./Fn. E4-PM13 E4-ESD E4-PM1 E4-PM3 E4-RQ E4-EGS Polar-PM13 Polar-ESD Polar-PM1 Polar-PM3 Polar-RQ Polar-EGSD Real-DE
Attractive Sector Median σ 2.589e-14 3.386e-13 1.467e-08 2.062e-05 2.641e-33 1.588e-32 1.176e-33 7.216e-33 0.456 0.344 0.0102 0.05168 2.824e-16 3.436e-14 4.38e-10 6.177e-04 4.584e-34 1.327e-33 1.393e-34 3.08e-34 0.02154 0.02112 2.492e-03 2.413e-03 1.024e-08 1.118e-06
Step Ellipsoidal Median σ 0 0 0 0.02647 0 0 0 0 0.09125 0.06045 0.14 0.09471 0 0 0 0 0 0 0 0 0.005246 0.01144 0.01211 0.04367 0 6.449e-03
in the benchmark are: the Rastrigin, the Weierstrass, the Schaffers F7, the moderately ill-conditioned Schaffers F7 and the Composite Griewank-Rosenbrock functions. 5) Multi-modal functions with weak global structure (MWeak). These functions have multiple local optima, but lack of a global pattern that guides the optimization algorithm toward the global optimum. Examples included in the benchmark are the Schwefel, the Gallagher’s Gaussian 101-me Peaks, the Gallagher’s Gaussian 21hi Peaks, the Katsuura and the Lunacek bi-Rastrigin functions. Due to the isomorphism between R3 and the imaginary part of the quaternion, we evaluated the algorithms with 3dimensional functions. In addition, the result was compared with the real-valued DE algorithm. The methodology of the experimental analysis was as follows: For each fitness function, we run 20 experiments using each algorithm and saved the fitness value of the best element of the population at each generation. Then we computed the the mode and standard deviation of the 20 experiments for each algorithm; these values are shown at Tables II to VI for each function of the benchmark. From this data, we selected the fitness value of the best individual in the last generation and applied a statistical analysis for testing four hypotheses: which scheme produces the best solution for all the problems on the benchmark,
Rosenbrock Median σ 9.585 9.143 0.207 2.12 5093.751 3501.937 3493.714 2989.844 1168.673 1652.548 1.454 2.233 186.662 3693.806 17.138 1043.506 8262.610 4561.445 11871.863 5678.366 7298.722 3615.84 2.394 2.472 0.04961 0.221
Rotated Rosenbrock Median σ 30.014 22.039 0.252 69.783 6346.729 11561.102 5816.107 7063.83 3529.338 4771.858 1.464 3.885 23.611 47.022 16.594 381.642 8573.18 15835.696 28794.927 13120.706 11763.882 12154.138 1.284 2.795 0.02598 0.895
which scheme is the best for each group of functions, which mutation function produces better results, and which one is the best initialization method. The statistical analysis consisted of applying the Friedman test with the corresponding Nemenyi post-hoc test with an alpha of 0.05 [23]. Figure 3 shows the resulting critical difference diagrams for the first two hypotheses. For the first test, the best results are produced by PM3 and PM1 mutation functions with any type of initialization method; in addition, the fitness of the real-valued DE algorithm is not is not statistically significantly different from the one obtained by these quaternion algorithms. However, if we analyze the results by the type of functions (second test), we can see that the best quaternion algorithms (PM1, PM3) outperform the real-valued versions for U-High functions as well as for Madequate functions. The poor performance of PM1 and PM3 in U-Low problems, is due to their poor performance on the Rosenbrock functions; thus, it can be concluded that the search process of these algorithms is not capable of follow a long path with D-1 changes in the direction. Moreover, the poor performance of PM1 and PM3 in M-Weak problems is due to their poor performance in all the functions, but in the Katsuura one. This can be explained by the fact that this function has more than 10D global optima; thus, the quaternion algorithms are capable of finding the optimal values in cases in which they
7
TABLE IV M EAN AND STANDARD DEVIATION OF THE FITNESS VALUES OBTAINED BY EACH ALGORITHM AT GENERATION 100 FOR FUNCTIONS WITH HIGH CONDITIONING AND UNIMODAL . Alg./Fn. E4-PM13 E4-ESD E4-PM1 E4-PM3 E4-RQ E4-EGS Polar-PM13 Polar-ESD Polar-PM1 Polar-PM3 Polar-RQ Polar-EGSD Real-DE
Rotated Ellipsoidal Median σ 0.08855 0.591 9.396e-08 3.621 8.559e-15 1.126e-18 8.564e-15 7.143e-18 36.522 73.296 234.085 775.839 0.01936 0.145 2.671e-07 0.03088 8.559e-15 3.041e-19 8.569e-15 1.351e-17 1.162 1.885 17.759 157.307 4.039e-03 0.03186
Discus Median σ 1.114e-14 1.922e-11 4.538e-05 1.565e-03 8.55e-15 5.881e-19 8.551e-15 8.949e-18 0.817 7.461 40.83 109.567 8.696e-15 5.127e-10 9.545e-08 0.162 8.55e-15 3.128e-20 8.551e-15 6.249e-18 0.06328 0.175 4.324 11.516 3.701e-03 0.152
Bent Cigar Median σ 17.299 585.498 0.144 1.585 8.604e-26 1.258e-04 1.969e-29 5.859e-26 380.022 1550.696 641.17 3300.459 3.327 30.278 0.03747 0.737 2.561e-27 4.184e-04 1.307e-32 4.116e-30 22.479 80.123 221.014 473.756 0.547 0.755
Sharp Ridge Median σ 0.735 5.99 8.17e-05 0.222 4.166e-15 4.21e-13 1.76e-15 8.246e-15 3.133 3.134 10.237 5.739 0.129 1.442 1e-05 2.862e-03 3.229e-14 3.459e-03 1.956e-16 1.917e-15 0.54 3.713 4.178 4.281 0.01167 0.02375
Different Powers Median σ 6.325e-09 3.409e-07 4.184e-08 0.04493 5.041e-22 4.826e-15 1.064e-18 9.807e-17 0.361 0.218 9.492e-03 6.838e-03 3.227e-11 2.481e-06 1.738e-08 2.124e-06 8.112e-22 2.528e-13 9.193e-18 2.009e-16 8.359e-04 8.876e-03 1.498e-03 3.911e-03 2.147e-06 7.471e-06
TABLE V M EAN AND STANDARD DEVIATION OF THE FITNESS VALUES OBTAINED BY EACH ALGORITHM AT GENERATION 100 FOR MULTI - MODAL FUNCTIONS WITH ADEQUATE GLOBAL STRUCTURE . Rastrigin
Weierstrass
Schaffers F7
Alg./Fn. E4-PM13 E4-ESD E4-PM1 E4-PM3 E4-RQ E4-EGS Polar-PM13 Polar-ESD Polar-PM1 Polar-PM3 Polar-RQ Polar-EGSD Real-DE
Median 1.581 4.669 7.274e-20 5.275e-20 3.671 6.358 0.416 0.632 7.065e-20 5.455e-20 1.694 0.262 2.332
σ 1.334 2.585 1.664 0.67 1.848 2.923 1.374 1.577 1.834e-20 1.698e-20 0.98 1.828 1.263
Median 1.17 1.824 3.832e-34 3.795e-34 1.863 2.446 1.372e-06 0.104 3.883e-34 3.821e-34 2.027 0.853 1.894
σ 0.647 1.767 1.847 6.939e-36 1.709 1.247 0.905 0.887 2.579e-35 6.913e-36 2.14 1.754 1.195
Median 3.217e-07 1.187e-03 1.249e-18 1.013e-18 0.406 0.201 2.488e-09 2.005e-03 1.304e-19 1.019e-19 0.05822 0.06018 7.289e-03
σ 5.702e-04 0.04699 3.638e-18 1.713e-17 0.282 0.177 3.23e-06 0.03086 2.693e-18 1.921e-19 0.143 0.07052 5.207e-03
Schaffers F7 Mod. Ill-conditioned Median σ 0.01751 0.182 0.03182 0.05982 2.513e-16 2.948e-13 8.422e-17 2.022e-15 1.164 0.576 1.428 1.002 9e-03 0.03529 0.02152 0.08615 1.603e-16 1.423e-05 1.677e-18 2.977e-17 0.139 0.134 0.447 0.391 0.02232 0.175
Composite Griewank-Rosenbrock Median σ 0.249 0.137 0.489 0.417 0.57 0.407 0.426 0.317 0.245 0.311 0.541 0.389 0.23 0.271 0.293 0.375 0.938 0.421 0.481 0.173 0.416 0.261 0.372 0.206 0.237 0.281
TABLE VI M EAN AND STANDARD DEVIATION OF THE FITNESS VALUES OBTAINED BY EACH ALGORITHM AT GENERATION 100 FOR MULTI - MODAL FUNCTIONS WITH WEAK GLOBAL STRUCTURE . Alg./Fn. E4-PM13 E4-ESD E4-PM1 E4-PM3 E4-RQ E4-EGS Polar-PM13 Polar-ESD Polar-PM1 Polar-PM3 Polar-RQ Polar-EGSD Real-DE
Schwefel Median σ 2.703 163.44 1.632 5125.971 20148.855 7464.203 14113.66 7049.684 8893.436 4712.046 1.597 0.508 238.573 4338.006 3.679 3303.669 34683.961 2915.314 39096.476 3261.3 30097.834 1288.747 1.801 0.421 0.463 0.373
Gallagher’s Gaussian 101 Median σ 1.415 1.142 0.06854 0.581 4.926 5.48 4.263 2.744 2.616 1.108 0.931 0.515 2.948 2.558 2.2 2.061 7.741 5.041 7.617 5.636 3.551 5.604 0.06162 0.93 0.718 0.727
Gallagher’s Gaussian 21 Median σ 2.04 5.503 0.07732 2.117 9.888 14.065 7.771 7.398 3.998 15.132 0.371 1.872 1.44 8.072 2.145 5.196 13.115 12.147 20.354 27.839 25.254 27.479 0.41 0.768 0.37 2.13
Katsuura Median σ 2.173 0.664 2.488 1.265 3.139 1.517 0.996 0.908 1.757 0.643 1.812 0.652 1.403 0.518 1.247 0.611 1.852e-08 1.417 6.284e-14 0.595 1.279 0.581 2.105 0.74 2.017 0.466
Lunacek bi-Rastrigin Median σ 7.293 1.725 5.863 3.564 16.978 3.85 12.924 3.971 5.954 2.914 9.902 3.579 8.109 3.326 8.241 3.943 18.742 1.252 17.892 1.235 10.348 3.343 6.655 2.478 6.237 2.421
8
All functions from the benchmark
S eparable functions
Low or moderate conditioning functions
13 12 11 10 9 8 7 6 5 4 3 2 1
13 12 11 10 9 8 7 6 5 4 3 2 1
13 12 11 10 9 8 7 6 5 4 3 2 1
5.1844 E 4-PM3 5.3979 Polar-PM3 5.5083 Polar-PM1 5.6188 E4-PM1 5.9979 Real-DE 6.2125 Polar-PM13 6.5667 Polar-ESD
E4-RQ 10.0292 Polar-RQ 9.4656 E4-EGSD 9.3833 Polar-EGSD 8.0573 E4-PM13 6.8385 E4-ESD 6.7396
E4-RQ 11.01 E4-EGSD 10.26 Polar-RQ 9.54 Polar-EGSD 8.65 E4-ESD 7.815 E4-PM13 7.41
Multi-modal with adequate global structure functions
High conditioning and unimodal functions
13 12 11 10 9 8 7 6 5 4 3 2 1
13 12 11 10 9 8 7 6 5 4 3 2 1
2.345 Polar-PM1 2.55 Polar-PM3 2.61 E4-PM1 2.93 E4-PM3 6.56 Polar-ESD 6.78 E4-ESD 6.93 Polar-PM13
E4-EGSD 12.31 E4-RQ 11.83 Polar-EGSD 11.07 Polar-RQ 9.64 Real-DE 7.86 E4-PM13 7.585
4.115 Polar-PM1 4.12 E4-PM1 4.205 E4-PM3 4.4 Polar-PM3 5.85 Real-DE 6.345 Polar-PM13 7.28 Polar-ESD
E4-EGSD 10.83 E4-RQ 9.95 Polar-RQ 8.78 Polar-EGSD 8.7 E4-ESD 8.7 Real-DE 8.02
3.375 Polar-PM3 3.635 E4-PM3 4.175 Polar-PM1 4.785 E4-PM1 5.7 Polar-PM13 7.14 E4-PM13 7.21 Polar-ESD
Polar-RQ 10.545 E4-RQ 10.29 Polar-PM1 8.445 Polar-PM3 8.315 E4-PM1 7.645 E4-PM3 7.495
4.055 R eal-DE 4.72 E4-ESD 5.26 E4-PM13 5.355 Polar-ESD 5.57 Polar-PM13 6.145 Polar-EGSD 7.16 E4-EGSD
Multimodal with weak global structure functions 13 12 11 10 9 8 7 6 5 4 3 2 1
Polar-PM1 9.86 E4-PM1 9.74 Polar-PM3 9.71 Polar-RQ 9.13 E4-PM3 8.46 E4-RQ 6.75
3.715 R eal-DE 4.82 Polar-EGSD 5.18 E4-ESD 5.19 E4-EGSD 6.045 Polar-ESD 6.18 E4-PM13 6.22 Polar-PM13
Fig. 3. Critical difference diagram showing the ranking of the QDE schemes based on fitness values, both across all functions and within each function group of the BBOB benchmark.
TABLE VII QDE SCHEMES THAT OBTAINED A FITNESS VALUE WITH NOT STATISTICALLY SIGNIFICANTLY DIFFERENCE FROM THE BEST RANKED ALGORITHM FOR EACH GROUP OF FUNCTIONS .
QDE E4-PM1 Polar-PM1 E4-PM3 Polar-PM3 E4-PM13 Polar-PM13 E4-ESD Polar-ESD E4-EGSD Polar-EGSD Real-valued
Sep. X X X X
U-Low
U-High X X X X
M-adequate X X X X
X X X X
X
M-weak
X X X X
X
Ranking of Quaternion Mutation Functions 6
RQ 4.6578 EGSD 4.2219 ESD 3.2984
5
4
3
2
1
2.724 PM3 2.8833 PM1 3.2146 PM13
Fig. 4. Critical difference diagram showing the ranking of the QDE mutation schemes based on fitness values.
are numerous. Another interesting finding is that the best quaternion algorithms (PM1, PM3) were outperformed by other quaternion algorithms for U-Low and M-Weak functions. In particular, the E4-ESD algorithm is capable of solving both type of problems. These results are summarized in Table VII. The third test consisted in comparing the quaternion-valued mutation functions; according to the statistical analysis, the best results are obtained by PM3 and PM1 functions, see Figure 4. For the last hypothesis, there is not statistically significantly difference between the results produced by the initialization methods.
Finally, we analyzed the convergence of the algorithms. To do this, we selected the last generation at which the fitness value ceased to change. Then, a statistical analysis was performed using the Friedman test, followed by the Nemenyi post-hoc test with a significance level of α = 0.05 [23]. Figure 5 presents the resulting critical difference diagram based on convergence data from all fitness functions. The diagram indicates that most quaternion-valued algorithms converge faster than the real-valued version. Additionally, the realvalued and the quaternion-valued variants with the highest accuracy–namely, Polar-PM3 and PM1 with any initialization method–did not show a statistically significant difference in the number of generations required for convergence. However, when analyzing the convergence of the top-ranked algorithms within each function group, we observed that the real-valued algorithm converged faster on U-High and MAdequate functions, which are the same functions where it performed worse in terms of fitness. In contrast, for Separable, ULow, and M-Weak functions, several QDE algorithms achieved faster convergence than the real-valued version, see Figure 5. Notably, Polar-PM1 and Polar-PM3 demonstrated both strong convergence and consistent fitness performance on Separable, U-High, and M-Adequate functions. For U-Low functions, PM13 with any initialization method performed best, while for M-Weak functions, the EGSD variant with any initialization method showed superior convergence. IV. D ISCUSSION In this section, we discuss some interesting findings when evaluating the performance of the QDE algorithm. Firstly, because of the No Free Lunch theorems [24], [25], no optimization algorithm is universally superior across all problems. Consequently, quaternion-valued algorithms offer advantages for specific groups of functions. Based on our statistical analysis, QDE algorithms using mutation strategies PM1 and PM3 with polar initialization are the most
9
All functions from the benchmark 13 12 11 10 9 8 7 6 5 4 3 2 1
Separable functions 4.921 E 4-R Q 4.934 Polar-RQ 5.118 E4-PM13 5.206 E4-EGSD 5.837 Polar-PM13 6.101 Polar-EGSD 7.758 Polar-PM1
Polar-ESD 9.225 E4-ESD 8.744 E4-PM3 8.585 Polar-PM3 8.35 E4-PM1 8.142 Real-DE 8.079
High conditioning and unimodal functions 5
E4-PM1 3.4 E4-PM3 3.255
4
3
2
1
5
4
2
Low or moderate conditioning functions
1 5
2.635 Polar-PM1 2.685 Polar-PM3 2.98 E4-PM3
Real-DE 3.415 E4-PM1 3.285
Multi-modal with adequate global structure functions 5
2.53 Real-DE 2.89 Polar-PM1 2.925 Polar-PM3
3
4
E4-PM3 3.46 E4-PM1 3.37
3
2
1
3
2
1
2.1 Polar-PM13
Polar-ESD 3.7875 E4-ESD 3.675
2.1875 E4-PM13 3.25 Real-DE
Multimodal with weak global structure functions 4
1.875 Real-DE 3.01 Polar-PM1 3.285 Polar-PM3
4
E4-ESD 3.125 Real-DE 2.84
3
2
1
1.95 E4-EGSD 2.085 Polar-EGSD
Fig. 5. Critical difference diagram showing the ranking of the QDE schemes based on convergence data, both across all functions and within each function group of the BBOB benchmark
effective when addressing U-High or M-adequate problems. These configurations not only achieve better optimization values but also demonstrate faster convergence. For separable functions, the real-valued algorithm achieves similar accuracy to these quaternion-based schemes; however, the quaternion algorithms are advantageous due to their faster convergence. When dealing with U-Low functions, PM13 with any initialization method yields the best performance in both fitness and convergence speed. In contrast, for M-Weak problems, the EGSD variant and the real-valued algorithm reach comparable accuracy, but the quaternion-based approach converges more quickly. Secondly, the faster convergence of quaternion algorithms could be explained by two factors: 1) the definition of the Hamilton product involves intercomponent products taking advantage of intercomponent relationships between variables, as shown in Equation 10, and 2) the geometric interpretation of Hamilton product causes a most suitable navigation on the search space. This also explains the superior performance of the Polar mutation algorithms, since multiplication of quaternions represents a geometrical rotation. This does the search of the optimum value more efficient than when applying arithmetic operations like the E4 Mutation method does. Next, we summarize some interesting behaviors when comparing similar algorithms: The difference between Simple Mutation and PM1 algorithms is that the first one uses a random quaternion, while the second constructs a quaternion from the population, this change causes a directed search of the PM1 algorithm, which outperforms the Simple Mutation algorithm. • PM1 and PM3 algorithms differ only in the use of a third individual for applying the resulting transformation. It could be argued that this would cause a better performance of algorithms PM3, due to adding an extra trait for avoiding local minimums; however, both algorithms obtained similar performance. • PM1 and PM13 algorithms differ just in the sum of the •
quaternion q3 . This causes that PM1 performs better for Separable, U-High and M-adequate functions, while the PM13 algorithm performs better for U-Low functions. • The superior performance of EGSD algorithms for Mweak functions could be explained by the use of multiple Hamilton products, see Equation 19. Note that in this case we have three quaternions took from the population and a random quaternion. Even though M-Weak functions does not present clear pattern, the interchannel products and the use of multiple individuals might help to solve the problem. This information could be missleadinf for simple problems, which would explain the poor performance on other group of functions. V. C ONCLUSIONS AND F UTURE W ORKS This paper introduces a family of novel quaternion-valued differential evolution algorithms, which, for some specific group of optimization functions, are capable of outperforming the real-valued versions, and for the rest of the functions obtain similar results to the real-valued versions. In addition, quaternion algorithms have a faster convergence to the optimal value in the vast majority of the cases. Thus, the proposed algorithms PM1, PM3 (with any initialization method) should be the first option to test when dealing with an optimization problem. Finally, this paper opens the possibilities of endowing optimization algorithms with novel search methods by applying geometric properties embedded in the algebraic representation. Thus, future works should focus on extending the use of quaternion operators to other bioinspired optimization algorithms, extending the QDE algorithm to the N -dimensional domain using Clifford algebras, and exploring the use of different signatures of algebra for searching in non-Euclidean spaces. In terms of applications, future works should explore the use of QDE to relevant real-world optimization problems with 3 or 4 decision variables, such as robotic kinematics and dynamics, 3D pose estimation from images, mechanical and aerodynamic
10
design, to name a few. These domains naturally benefit from quaternion-based modeling due to their inherent spatial and rotational properties. ACKNOWLEDGMENTS G.Altamirano received a Postdoctoral Fellowship Estancias Posdoctorales por México from CONAHCYT. C ONFLICT OF INTEREST STATEMENT All authors declare that they have no conflicts of interest. R EFERENCES [1] P. Wang, Y. Zhou, Q. Luo, C. Han, Y. Niu, and M. Lei, “Complexvalued encoding metaheuristic optimization algorithm: A comprehensive survey,” Neurocomputing, vol. 407, pp. 313–342, 2020. [2] I. Fister, X.-S. Yang, J. Brest, and I. Fister Jr., “Modified firefly algorithm using quaternion representation,” Expert Systems with Applications, vol. 40, no. 18, pp. 7220–7230, 2013. [3] D. Saldias, L. Radavelli, C. Roesler, and D. Martins, “Kinematic synthesis of the passive human knee joint by differential evolution and quaternions algebra: A preliminary study,” in Proceedings of the 5th IEEE RAS/EMBS International Conference on Biomedical Robotics and Biomechatronics, ser. BioRob, Sao Paulo, Brazil, 2014, pp. 96–101. [4] I. Fister, J. Brest, I. Fister, and X.-S. Yang, “Modified bat algorithm with quaternion representation,” in Proceedings of the IEEE Congress on Evolutionary Computation (CEC), ser. CEC, Sendai, Japan, 2015, pp. 491–498. [5] T. T. Khuat and M. H. Le, “A genetic algorithm with multi-parent crossover using quaternion representation for numerical function optimization,” Applied Intelligence, vol. 46, no. 4, p. 810–826, 2017. [6] S. Song, X. Chen, Y. Zhang, Z. Tang, and Y. Todo, “Protein–ligand docking using differential evolution with an adaptive mechanism,” Knowledge-Based Systems, vol. 231, p. 107433, 2021. [7] N. Arana-Daniel, “Complex and hypercomplex-valued support vector machines: A survey,” Applied Sciences, vol. 9, no. 15, 2019. [8] T. Parcollet, M. Morchid, and G. Linàres, “A survey of quaternion neural networks,” Artificial Intelligence Review, vol. 53, no. 4, pp. 2957–2982, 2020. [9] G. Altamirano-Gomez and C. Gershenson, “Quaternion convolutional neural networks: Current advances and future directions,” Advances in Applied Clifford Algebras, vol. 34, no. 42, pp. 1–63, 2024. [10] Z. Zhou, Y. Chen, and Y. Zhou, “Simultaneously learning deep quaternion reconstruction and noise convolutional dictionary for color image denoising,” IEEE Transactions on Emerging Topics in Computational Intelligence, vol. 9, no. 2, pp. 1766–1779, 2025. [11] Y. Zhou, L. Jin, G. Ma, and X. Xu, “Quaternion capsule neural network with region attention for facial expression recognition in color images,” IEEE Transactions on Emerging Topics in Computational Intelligence, vol. 6, no. 4, pp. 893–912, 2022. [12] J. Liu, X. Liao, and J.-S. Dong, “A quaternion-valued neural network approach to nonsmooth nonconvex constrained optimization in quaternion domain,” IEEE Transactions on Emerging Topics in Computational Intelligence, vol. 8, no. 1, pp. 654–669, 2024. [13] Q. Yu, Z. Zhang, Z. Wang, H. Li, Z. Lei, and S. Gao, “Quaternion dendritic neuron model for multivariate financial time series prediction,” IEEE Transactions on Emerging Topics in Computational Intelligence, pp. 1–12, 2024. [14] G. Sfikas, G. Retsinas, A. Giotis, B. Gatos, and C. Nikou, “Keyword spotting with quaternionic ResNet: Application to spotting in Greek manuscripts,” in Document Analysis Systems: 15th IAPR International Workshop, ser. DAS, La Rochelle, France, 2022, pp. 382–396. [15] W. Hamilton, Lectures on quaternions: Containing a systematic statement of a new mathematical method. Dublin, England: Hodges and Smith, Whittaker & Co., MacMillan & Co., 1853. [16] ——, Elements of quaternions. London, England: Longmans, Green, & Co., 1866. [17] W. R. Hamilton, “On quaternions, or on a new system of imaginaries in algebra,” 2000, https://www.maths.tcd.ie/pub/HistMath/People/Hamilton/OnQuat/. [18] I. Kantor and A. Solodovnikov, Hypercomplex numbers. New York, NY: Springer-Verlag, 1989.
[19] J. Ward, Quaternions and Cayley numbers. Dordrecht, Netherlands: Kluwer Academic Publishers, 1997. [20] A. Hanson, Visualizing quaternions. San Francisco, CA: Elsevier Inc., 2006. [21] P. Du Val, Homographies, quaternions and rotations. London, England: Oxford University Press, 1964. [22] N. Hansen, A. Auger, R. Ros, S. Finck, and P. Pošı́k, “Comparing results of 31 algorithms from the black-box optimization benchmarking bbob-2009,” in Proceedings of the 12th Annual Conference Companion on Genetic and Evolutionary Computation, ser. GECCO ’10, Portland, Oregon, USA, 2010, pp. 1689–1696. [23] J. Demšar, “Statistical comparisons of classifiers over multiple data sets,” Journal of Machine Learning Research, vol. 7, pp. 1–30, 2006. [24] D. Wolpert and W. Macready, “No free lunch theorems for search,” The Santa Fe Institute, Santa Fe, NM, Tech. Rep. SFI-WP-95-02- 010, Feb. 1995. [25] D. H. Wolpert and W. G. Macready, “No free lunch theorems for optimization,” IEEE Transactions on Evolutionary Computation, vol. 1, no. 1, pp. 67–82, 1997.