{'fraction_non_alphanumeric': 0.035761702856653386, 'fraction_numerical': 0.012291746969891644, 'mean_word_length': 4.371520588911893, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Simple Yet Efficient Rank One Update for Covariance Matrix Adaptation

In this paper, we propose an efficient approximated rank one update for covariance matrix adaptation evolution strategy (CMA-ES). It makes use of two evolution paths as simple as that of CMA-ES, while avoiding the computational matrix decomposition. We analyze the algorithms' properties and behaviors. We experimentally study the proposed algorithm's performances. It generally outperforms or performs competitively to the Cholesky CMA-ES.


Introduction
Evolution strategies (ES) are one of the main streams of evolutionary algorithms [5] for continuous black-box optimization [13] [20]. Typically, ES evolves a Gaussian distribution [5]. It iteratively samples one or more new candidate solutions, and adapts the distribution parameters for producing more promising solutions. The covariance matrix adaptation evolution strategy (CMA-ES) is one of the most successful evolution strategies. It adapts a full covariance matrix, which achieves invariant under general linear transformation in the search space [7] [8].
The covariance matrix is adapted to increase the probability of reproducing the successful search directions in the previous generations. In practice, it is updated by the rank-µ update with maximum likelihood estimation of successful solutions and the rank one update with the evolution paths, which involves the historical search information. A main limitation of CMA-ES is that it has to conduct matrix decomposition to generate new solutions [10] [9] [6], which is computationally expensive. This prohibits its applications to problems with over hundreds of variables. In practice, the matrix decomposition can be done every O(n) generations to alleviate the computational overload, during which the outdated matrices are used to generate new solutions. This is undesirable as it does not make use of the current search information, and may result in some performance loss [23].
To alleviate the issue of computational overload, Cholesky CMA-ES directly adapts the Cholesky factor instead of the covariance matrix [12]. It considers a rank one perturbation of the covariance matrix, and obtains a rank one update for the Cholesky factor. This achieves complexity O(n 2 ). To use the evolution path in Cholesky CMA-ES, it has to maintain and update the inverse Cholesky factor accordingly [23], which is much complicated than the rank one update for the covariance matrix. It becomes even more complicated when restricting the Cholesky factor to be upper triangular [14].
In this paper, we propose to construct a new evolution path to approximate the inverse vector of the evolution path A −1 p. Then both evolution paths contributes an efficient rank one update for the mutation matrix A that satisfies AA T = C. It avoids maintaining the inverse Cholesky factor, and reduces the computational complexity in the update procedure of Cholesky CMA-ES to a half. The obtained algorithm is as simple as CMA-ES with pure rank one update [10].
We analysis the approximations in the proposed algorithm, and validate the approximations experimentally. We experimentally study the algorithm performances. It shows that the proposed methods outperform or performs competitively to Cholesky CMA-ES.

Background
In this section, we consider the black-box optimization in the continuous domain (1) In black-box optimization, we have no gradient information of f available, and the only information accessible is the objective value for a given point x ∈ R n .

Gaussian Mutations in Evolution Strategies
In a canonical ES, a new candidate solution is generated (at generation t) by where m t is the current distribution mean, determining the location of the mutation distribution over the search space, and σ t is the step size which determines the global variance of the underlying distribution. The mutation y ∈ R n is a n-dimensional random vector distributed according to the Gaussian mutation distribution y ∼ N (0, C t ).
The approaches that evolving arbitrary Gaussian mutation distributions introduce complete invariance with respect to translation and rotation in the search space [11].
To draw a mutation from the Gaussian distribution N (0, C t ), it is necessary to decompose the covariance matrix. In CMA-ES, the eigen-decomposition is conducted as C t = BD 2 B T , where B is orthogonal matrix with eigenvectors as columns, and D is a diagonal matrix with the square root of eigenvalues as the entries. A new solution is generated by An alternative is to use the Cholesky decomposition [6]. The computational complexity of the matrix decomposition is O(n 3 ) in general, which dominates the computational complexity of CMA-ES.

The Parameter Dependencies
Evolving a full covariance matrix renders the invariance under generalized linear transformations in the search space [7] [8]. The covariance matrix or its factorizations actually present a transformation in the search space, which transforms the objective function to be more separable and well-scaled [7]. The transformed objective function can be solved more efficiently. Hence, it is of crucial importance to learn the local transformation.

Covariance Matrix Adaptation
CMA-ES evolves a full covariance matrix to represent the parameter dependencies. It updates the covariance matrix by the maximum likelihood estimation of the successful mutations. The original CMA-ES [10] with pure rank one update is The evolution path p accumulates the movements of the distribution mean in consecutive generations and represents one of the most promising search directions [15].
To make better use of the current population, it also uses the search directions of the selected solutions. Thus, the update for covariance matrix is given by where c 1 , c µ are constant changing rates, and y i:λ denotes the mutation corresponding to the i-th best solution of the current population. The rank-µ update increases the variance of successful search directions by a weighted maximum likelihood estimation, and increases the probability of reproducing such directions.

Cholesky Covariance Matrix Adaptation
To avoid the costly matrix decomposition, Cholesky CMA-ES directly adapts the Cholesky factor of the covariance matrix C t = A t A T t . From the rank one update for the covariance matrix, it obtains the rank one Cholesky update [12].
To make use of the evolution path, it has to additionally maintain and adapt the inverse Cholesky factor A inv t to avoid the expensive computation of A −1 t p t+1 . This makes the adaptation scheme complicated [23]. The inverse vector of the evolution path is given by where A inv t approximates A −1 t . Then both matrices are updated by This update scheme reduces the internal complexity to O(n 2 ). However, it is much more complicated than the rank one update of CMA-ES (6).

Proposed Algorithm
We propose an efficient update scheme for the matrix A that satisfies C ≈ AA T , where A can be referred to mutation matrix. The algorithm maintains and evolves the following parameters at generation t: • m t : the distribution mean which determines the location of the sampling region.
• A t : the non-singular mutation matrix, determining the shape of the sampling distribution.
• σ t : the step size which determines the global variance of the sampling region.
• p t , v t : the evolution paths used to update the mutation matrix.
• s t : the evolution path used to update the step size.
For convenience, we denote the evolution paths as the p-path, the v-path, and the s-path, respectively. At the beginning, these distribution parameters are initialized to be m t ∈ R n , σ > 0, A 0 = I, and the evolution paths 0.

Main Procedures
In the following, we present the major procedures of the algorithm.

Sampling New Solutions
At each generation after initialization, a number of λ solutions are generated from the current distribution. A new solution x i , i = 1, · · · , λ is obtained by

Evaluation and Selection
The objective values of the newly sampled solutions are computed to evaluate their qualities. Then they are sorted according to the objective values where x i:λ denotes the i-th best solution of the current population. For truncation selection, the best µ ≤ λ 2 ones are selected to update the distribution parameters. The distribution mean is updated as where w i can be set to w i = 1/µ for i = 1, . . . , µ. A better choice is given by It is designed w 1 ≥ · · · ≥ w µ to favor the top ranked solutions, and

Evolution Paths
The matrix A is updated by evolution paths. For convenience, we denote where z i:λ corresponds to the i-th best solution x i:λ . Obviously, we have y w = A t z w . We construct and update the evolution paths as The involved parameters are given in the following.
• The factor µ eff = 1/ µ i=1 w 2 i normalizes the update direction such that under random selection, we have • The c(2 − c) is determined by This indicates that, given under the random selection. The similar condition holds for the v-path that given v t ∼ N (0, I), we have v t+1 ∼ N (0, I).
• The constant changing rate c determines the number of generations involved in the cumulation of the evolution paths. Generally, the changing rate c is set to c ∝ 1/n such that c −1 ∝ n, which is the same order to the number of free parameters of p and v.
The p-path accumulates average successful search directions in consecutive generations. Averaging over generations, opposite search directions are canceled out and consistent search directions are accumulated. Hence, the p-path represents the most promising search direction [15].
Further, for c = 1, no cumulations take place, and the equations (15) and (16) gives that v t+1 = A −1 t p t+1 . For c ∈ (0, 1), the cumulations take place in both evolution paths, and the v-path accumulates z w = A −1 t y w at each generation. Both p 0 and v 0 are initialized to zero, and the deviations of the v-path from the p-path in previous generations gradually damped by the factor (1 − c). Therefore, for c ∈ (0, 1), we have the following approximation in general

Adaptation of the Mutation Matrix
We use the v-path to replace the inverse vector u t+1 , and update the mutation matrix as where c 1 is constant changing rate for the mutation matrix. Note that p t+1 v T t+1 is a n × n matrix with rank one. We denote this update scheme by mutation matrix adaptation evolution strategy (MMA-ES).
This rank one update can be obtained by analyzing the coefficients of the rank one update of the efficient Cholesky CMA-ES (eCMA-ES) [16]. In eCMA-ES, it updates the matrix A by two evolution paths as Substituting the approximation v t+1 ≈ A −1 t p t+1 , it is easy to verify that It is the rank one update of the covariance matrix in CMA-ES.
The changing rate c 1 is set to c 1 = 2/(n + √ 2) 2 in the Cholesky CMA-ES [23]. It gives that c 1 1 for even moderate n, which gives the approximation Further, as the v-path is normally distributed v ∼ N (0, I) under random selection, its norm is distributed according to the Chi distribution v ∼ χ(n) and E v = √ n. Hence, we generally have c1 (1−c1) v 2 1 for v ∼ N (0, I). Denote the coefficient by b, we obtain the following approximation The approximations are due to √ 1 + x ≈ 1 + x 2 for |x| 1. From the above approximations, we obtain the approximated rank one update for the matrix A given by equation (20). The matrix A satisfies C ≈ AA T . Therefore, MMA-ES is a great simplification of the Cholesky CMA-ES. Fig. 2 shows the values of b, c1 2 , and their differences on different dimensions. We use the expectation of v in the computation of b. It is obvious that b and c1 2 are quite close, and their differences are about several orders of magnitude smaller than b and c1 2 . Hence, there would be only subtle performance difference between these update schemes.

Cumulative Step Size Adaptation
Apart from the distribution mean and mutation matrix, the step size is adapted independently following the cumulative step size adaptation [10] [1]. Another s-path is constructed to exploit the correlations between consecutive search directions over generations Then the step size is adapted by comparing the length of the evolution path with its expectation where d σ ≥ 1 is the damping parameter. The E[ N (0, I) ] is the expected length of random vectors drawn from the standard Gaussian distribution. If consecutive directions z w are correlated, the s-path would be very long, and the step size tends to increase. If they are anti-related, the s-path would be rather short, and the step size tends to decease. The s-path is an alternative to the evolution path in CMA-ES [10] that accumulates successive search directions without scaling in the local coordinate system expanded by the eigenvectors. In the proposed algorithm, it is computationally expensive to obtain the orthogonal matrix, and the s-path is used as an alternative.

Algorithm Framework
The proposed algorithm is presented in the following Algorithm 1. The parameters are initialized in line 2. At each generation, a number of λ solutions are generated in line 5-7. These new solutions are evaluated and sorted in line 9. The distribution mean and search directions are computed in line 10. Two evolution paths are constructed in line 11-12, and the matrix A is updated in line 13. Then the step size is adapted in line 14-15.
The update equation for the mutation matrix A involves only the outer product of the evolution paths. Hence, the computational complexity is of O(n 2 ). As the sampling procedure involves matrix-vector multi-plication, it cannot be further reduced.

Parameter Settings
The involved algorithm parameters are borrowed from [23], and listed in Table 1. The changing rate c, c 1 are designed to be inverse proportional to the number of free parameters in the evolution paths and the matrix. Therefore, the number of generations involved in the cumulations of A is 7: x i = m t + σ t y i 8:

Analysis on the Evolution Path
We analyze the orthogonality of consecutive z-vectors of the v-path under some reasonable assumptions. For convenience, we denote the √ µ eff z w at generation t by z t . As previously demonstrated, the v-path is distributed according to the standard Gaussian distribution under random selection. Hence, we make the assumptions that v t and z t are unchanged by selection over generation.
That is, only the directions of v t and z t change while their lengths are kept constant. By substituting (16) into v t+1 ≈ v t , we can get v T t z t+1 ≈ 0. Given (c 1 + c µ ) 1, we have A t ≈ A t+1 . By further assuming v T t−1 z t+1 ≈ 0, we can derive that This means that consecutive mutation steps z t and z t+1 tends to be orthogonal. This can be equivalently written as where C t = A t A T t . This indicates that successive search directions are conjugate in terms of the covariance matrix. As the covariance matrix generally converges to the inverse Hessian matrix on quadratic objective functions (up to a scalar factor) [10] [3], this is desirable for efficient search.

Experimental Studies
We investigate the algorithm behaviors and performances on a set of commonly used test problems defined in table 2. We focus on the convergence properties of the algorithm. The initial distribution mean m 0 is uniformly randomly sampled in [−10, 10] n , and the step size σ 0 is initialized to be 20/3. The algorithm runs 21 times independently on each test problem. In the comparison of the experimental results, we conduct Wilcoxon rank sum test to assess the difference. When we state that two results are different, we have assert its statistical significance at the 5% significance level.

Experiments on the Approximations
We experimentally validate the effectiveness of the approximation of the v-path to the inverse vector of the p-path, and the approximation of coefficients.

Approximation to the Inverse Vector
We first conduct experiments to validate the approximation of the v-path to the inverse vector u. As the v-path is distributed according to the standard Gaussian distribution under the stationarity condition, we care only the direction of v. At each generation, we compute the v t+1 using (16) and the inverse vector u t+1 = A −1 t p t+1 . The similarity between u t+1 and v t+1 is measured by It indicates that v is close to u if α ≈ 0. The experiments are conducted on n = 32. Fig. 3 presents the experimental results on some test problems. It shows that α is of the order 10 −3 to 10 −2 in the whole optimization procedure on the test problems, including the f Ros function which has a parabolic shaped valley. Consequently, the v-path approximates well to the inverse vector u.
We investigate the α on different dimensions. On each dimension, it is averaged over 21 independent runs. The experimental results are shown in Fig. 4. Clearly, α is of the order 10 −3 to 10 −2 on all dimensions of the test problems, including on f Ros which has a parabolic shaped valley. Hence, the v-path can be used as a good approximation to the inverse vector.

Approximation of the Coefficient
We then investigate the difference between |b − c1 2 | in MMA-ES with b given by (??). At each generation, we compute |b − c1 2 |. Fig. 5 presents the experimental results on n = 32. It clearly shows that the difference |b − c1 2 | is about the order of 10 −5 in the whole optimization procedure, indicating that the rank one update (20) is a good approximation to the Cholesky update.
(a) f Cig  We further investigate the difference |b − c1 2 | on different dimensions, which is averaged over the whole optimization procedure of 21 runs on each dimension. The experimental results in Fig. 6 show that |b − c1 2 | descends from 10 −2 to 10 −8 as the dimension increases from of n = 2 to n = 256. Hence, the MMA-ES provides a good approximation to the Cholesky update.
(a) f Cig

Algorithm Behaviors
We investigate the algorithm behaviors in this subsection, including the evolution of the eigenvalues of the covariance matrix, and the performance invariance under rotation in the search space.

Evolution of the Eigenvalues
We investigate the evolution of the eigenvalues of C = AA T in the optimization procedure. At each generation, we decompose to C t = BD 2 B T , where D = diag(d 1 , . . . , d n ) with d i denoting the square root of the i-th eigenvalue of C t . d i represents the scale of the variables along each eigenvector. Fig. 7 shows the experiential result on f Ell with n = 10. The algorithm gradually learns the relative scales along the eigenvectors in the optimization procedure, and once the mutation matrix learns the relative scales, the objective value and step size decrease rapidly.

Invariance Under Rotation in the Search Space
The main advantage of CMA-ES is its invariance properties under general linear transformation in the search space, especially the rotation. We investigate the invariance under the rotation of the search space. We generate a orthogonal matrix R via the Gram-Schmidt process, and investigate the performance of the algorithm on f (Rx). Fig. 8 shows the median run out of 21 independent runs on f Ell and f Ros with both n = 32 and n = 64. It clearly shows that the algorithm performs uniformly on the f (x) and f (Rx). The slight differences are mainly due to the initialization of the distribution mean m 0 . Hence, the algorithm performance remains invariant under rotation in the search space.

Algorithm Performances Comparison
we compare the algorithms performances on the test problems. We present the median run in terms of the number of function evaluations of each algorithm. As the f Ros function has a local minimum [22] and the algorithms may fail, we consider that the failed runs cost more objective function evaluations than any successful runs. Fig. 9 shows the experimental results. The most important observation is that MMA-ES performs quite well on all the test problems. It outperforms Cholesky CMA-ES on most of the test problems with ill-conditioned Hessian matrix, and the f pR with a parabolic ridge. This validates the effectiveness of MMA-ES update. On the f Sp , MMA-ES performs similarly to Cholesky CMA-ES. On the f Ros , MMA-ES costs slightly more function evaluations than that of Cholesky CMA-ES. This may due to the parabolic shaped valley in the fitness landscape of f Ros , which requires rapid adaptation of the evolution path and the v-path may lay behind. MMA-ES performs quite similar on the problems, while the update scheme is much simpler.

Scaling on Dimension
We investigate the algorithm scalability on dimension of the search space. The algorithm scalability is measured by the average number of objective function evaluations (FE) on dimension. Because the f Ros is multi-modal [22], we use the success performanceŜP = E s /p s , where E s is the average function evaluations of successful runs, and p s is the success ratio of all independent runs [2] [11]. The experimental results in Fig. 10 and Fig. 11 show that the proposed method scales similarly to the Cholesky CMA-ES, while presenting some advantages on ill-conditioned problems. This indicates that, with increasing problem dimension, the proposed simplified update scheme works well.
To make the difference more clear, we investigate the relative performance of the proposed algorithms, with the Cholesky CMA-ES as baseline. The relative performance is measured by where F E 1 denotes the average number of function evaluations of the Cholesky CMA-ES, and F E 2 denotes that of the MMA-ES and the eCMA-ES. It measures how much times of function evaluations are required for the compared algorithms to reach the predefined accuracy. If F E2 F E1 < 1, the compared algorithm costs less function evaluations than Cholesky CMA-ES. The experimental results show that F E2 F E1 < 1 on most of the test problems and dimensions. Actually, Cholesky CMA-ES only outperforms MMA-ES and eCMA-ES on some dimensions on f Ros and f Sch .

Conclusion
In this paper, we have proposed an efficient rank one update for Cholesky CMA-ES by using a new v-path. It is as simple as the rank one update of CMA-ES, while avoiding the computational matrix decomposition. The v-path accumulates z-vectors drawn from the standard Gaussian distribution, and approximates the inverse vector of the p-path.
We have analyzed the orthogonal properties between consecutive z-vectors. We have conducted experiments to verify that the effectiveness of the approximations, and investigated the algorithm behaviors. Then we have showed that it outperforms or performs competitively to Cholesky CMA-ES. We have investigated the algorithm scalability on dimension.
The proposed method can be directly extended to the limited memory variant to handle large scale optimization problems following the idea of [18] [17]. The sparsity technique [19] can be used here without conducting matrix decomposition. The update (20) (especially for the case without cumulation in the evolution paths with c = 1) shows some similarity to the Hebbian learning rule [21]. The relation is to be further explored.
(a) f Ros Figure 11: The average number of function evaluations (FE) required to reach the predefined target, and the relative performance of two algorithms (II). The algorithms run 21 times independently on each dimension of the problems.