A covariance matrix adaptation evolution strategy in reproducing kernel Hilbert space

The covariance matrix adaptation evolution strategy (CMA-ES) is an efficient derivative-free optimization algorithm. It optimizes a black-box objective function over a well-defined parameter space in which feature functions are often defined manually. Therefore, the performance of those techniques strongly depends on the quality of the chosen features or the underlying parametric function space. Hence, enabling CMA-ES to optimize on a more complex and general function class has long been desired. In this paper, we consider modeling the input spaces in black-box optimization non-parametrically in reproducing kernel Hilbert spaces (RKHS). This modeling leads to a functional optimisation problem whose domain is a RKHS function space that enables optimisation in a very rich function class. We propose CMA-ES-RKHS, a generalized CMA-ES framework that is able to carry out black-box functional optimisation in RKHS. A search distribution on non-parametric function spaces, represented as a Gaussian process, is adapted by updating both its mean function and covariance operator. Adaptive and sparse representation of the mean function and the covariance operator can be retained for efficient computation in the updates and evaluations of CMA-ES-RKHS by resorting to sparsification. We will also show how to apply our new black-box framework to search for an optimum policy in reinforcement learning in which policies are represented as functions in a RKHS. CMA-ES-RKHS is evaluated on two functional optimization problems and two bench-marking reinforcement learning domains.


Introduction
The covariance matrix adaptation evolutionary strategy (CMA-ES) is a derivative-free method [12] which is a practical optimization tool for continuous optimization problems. It is a general optimization framework that possesses many appealing characteristics, e.g. derivative-free, covariant, off-the-shelf, scalable etc. It is especially useful on problems that are non-convex, non-separable, illconditioned, multi-modal, and with noisy evaluations. As the name indicates, CMA-ES belongs to Evolution Strategy (ES), hence it also works by operating three major steps: recombination, mutation and selection. In the context of robotics, CMA-ES has been widely used in many tasks: biped locomotion [43], whole-body locomotion optimization [8,9], swimming by [34], skill learning via reinforcement learning [14,15,31], inverse reinforcement learning [5,28], etc. Applying CMA-ES requires explicitly a finite-dimensional search space on which solution candidates live. In many domains, e.g. robotics, an optimization objective is often defined as a cost function of another parametric solution function. For instance, it might be an overall cost function depending on a robot controller, e.g. applications CMA-ES in robot skill learning [31] and policy search [15], or a loss function in the context of inverse optimal control [5,28], etc.. A robot controller is usually a parametric function of predefined feature maps, e.g. radial basis function (RBF) features. Though showing a lot of remarkable successes, such parametric approaches as well as their black-box solver, the parametric CMA-ES, could not model a very rich and flexible solution space, e.g. a complex behavior space. The performance of CMA-ES is affected not only by a suboptimal choice of features, but also the amount of used features. Therefore, integrating the ability of adaptive feature selection or non-parametric techniques into CMA-ES has been desired and not yet investigated, though this integration have been seen very often in other machine learning algorithms.
In this work, we propose CMA-ES-RKHS that enables functional optimization over a non-parametric solution space. Specifically, we assume that the solution (domain) space is a reproducing kernel Hilbert space (RKHS). Each solution candidate is a function in RKHS. Modeling the solution space this way, CMA-ES-RKHS is able to not only inherit full characteristics from CMA-ES, but enjoy other appealing properties of kernel methods as well. Firstly, CMA-ES-RKHS is able to optimize a functional objective whose domain consists of functions in a RKHS with an embedding kernel. That means the solution space does not need to depend on a set of predefined features that often has limited complexity, therefore our method can have richer representation power. Secondly, by modeling the solution space in RKHS, all update steps in CMA-ES-RKHS are handled analytically. We show that updated mean functionals, other intermediate terms, evolution path functionals or conjugate evolution path functionals are functions in the underlying RKHS. Moreover, the updated covariance is also an operator on the underlying RKHS. In other words, there is no additionally introduced computation when compared to the standard CMA-ES method. Thirdly, via sparsification techniques, a very complex search space can be represented compactly due to the reproducing property of the RKHS. However we can still achieve a solution of guaranteed quality which is at least comparable to that of the standard CMA-ES method. Though the use of sparsification might lead to representation with only a finite number of features, it allows adaptive feature selection and offers compact representation that reduce the computation cost.
By employing CMA-ES-RKHS we propose a non-parametric direct policy search technique in which the policy space in RL is an RKHS. As demonstrated by [2,19,36,38,40], the policy in RKHS is shown to be powerful and flexible in solving complex RL tasks. However their approaches use functional gradient ascent to update the policy functional that might have a problem with step-sizes or local optimality. It means the optimized policy does not fully exploit the flexible power of modeling in RKHS, which might consist of any complex policies including the global optimum one. As a global optimization method, our policy search via CMA-ES-RKHS is expected to search for such globally optimal policies. This paper is an extension from our work published [39].
Recently, there have been similar efforts in proposing new machine learning frameworks for functional optimization, such as functional regression by [16], representing motion trajectories in RKHS by [6,21], or finding geodesic shortest paths in physical systems by [17]. We demonstrate the proposed CMA-ES-RKHS and direct policy search via CMA-ES-RKHS in four experiments. The first two synthetic domains demonstrate the advantages of functional optimization in RKHS. The last two experiments are two RL domains, inverted pendulum and double-pendulum swing-up. The experiments will evaluate how CMA-ES-RKHS can do adaptive feature selection in a non-parametric way.
Our paper is organized as follows. The Background section will briefly discuss necessary preliminaries that will be used in the rest of the paper. Section 3 will start with a problem statement of functional optimization and then we will present our proposed framework, CMA-ES-RKHS. Section 4 will describe an application of CMA-ES-RKHS to a policy search problem in reinforcement learning. Section 5 will present experiment results on numerical and RL domains. Finally, we will conclude the paper with some remarks and future research directions.

Background
We briefly present a background of the Covariance Matrix Adaptation -Evolution Strategy (CMA-ES), its application for direct policy search reinforcement learning, and the recently introduced policy search in RKHS.

Covariance matrix adaptation-evolution strategy
CMA-ES is a global optimisation method introduced by [12]. It works by forming a parametric distribution over the solution space, e.g. the space of policy parameter in policy search, or the space of parameters of the loss function in inverse optimal control, etc. It iteratively samples a population of solution candidates from a parametrized search distribution. These candidates are then evaluated by a black-box function. Tuples of candidate-evaluation make up a dataset in order for CMA-ES to update the search distribution, i.e. its mean and its covariance matrix. Specifically, a cost function f ∶ ℜ n ↦ ℜ is parametrized by a parameter space ∈ ℜ n , f ( ) . The target is to find an optimum parameter * such that f ( * ) is minimum. It is common that a CMA-ES algorithm maintains a multivariate Gaussian distribution over the solution space as ∼ N( ; , ) , where is a n-dim mean vector and is a n × n covariance matrix. At each iteration k, it generates the kth population of offsprings from the kth distribution as i ∼ N( ; k , k ) , i = 1, ⋯ , , where k , k denote the mean vector and the covariance matrix at iteration k (after and got updates k times). Then, the offsprings are sorted ascendingly according to their evaluations f ( i ) . Only the first ( < ) best candidates are selected for use in updates of k and k . Another parameter is the global step-size ∈ ℜ that controls the convergence rate of the covariance matrix update. The parameter is defined as a global standard deviation. Hence, a full set of parameters in CMA-ES is { , , }.
In Algorithm 1, we give a full summary of the CMA-ES algorithm. Algorithm 1 starts with an initialisation of the parameters in steps 1 and 2. As discussed above, the parameters { , , } are updated iteratively as shown in steps from 4 to 13.
Step 4 is sampling from a normal distribution with mean and covariance . The evaluations via query to the black-box function f (⋅) are in step 5. The updated mean is a weighted sum of the best candidates as in step 7, in which the weights w i are set to 1∕ or to a better choice log( ∕2) − log(i) . The notation i∶ means the best candidate out of i , … , . The covariance matrix update in step 13 consists of three parts: (1) old information, (2) rank-1 update which computes the change of the mean over time as encoded in the evolution path c , and (3) rank-update which takes into account the good variations in the last population.
Step 10 devotes to the stepsize control update that constrains the expected changes of the distribution. Thus, this step is based on the conjugate evolution path . It aims to accelerate convergence to an optimum, and meanwhile prevents premature convergence. The other parameters: w is the variance effective selection mass, c 1 , c c , c are learning rates, and d is a damping factor for . The setting of these parameters is well studied and discussed in-depth in [10]. The termination can be set in various ways depending on the contracting covariance, or fitness functions (e.g. if the fitness functions of the population do not change for some iterations).
The updates of CMA-ES can alternatively be derived using the information-geometric concept of a natural gradient as shown in [11], which shares the same insight with the natural evolution strategies (NES) [44].
There are less efficient but practically fast techniques that can also adapt the covariance matrix: estimation of distribution algorithms (EDA) and the cross-entropy method (CEM). The major difference from CMA-ES is the choice of the reference mean value. EDA and CEM estimate the variance using the current population information, 1∶ , instead of exploiting old information as encoded in previous information of and c . Specifically, for the Gaussian search distribution, one can modify steps 9 and 13 in Algorithm 1 to receive new updates at iteration k for EDA and CEM [26,27]

as follows
The difference between EDA and CEM is: EDA updates C to the empirical covariance matrix, meanwhile CEM updates C to the unbiased empirical covariance matrix

Direct policy search in reinforcement learning
In this section, we describe the background of Markov decision process, reinforcement learning and policy search techniques.

Markov decision process and reinforcement learning
A Markov Decision Process (MDP) is a mathematical framework often used to formulate sequential decision making problems that is understood as an interaction process between an agent and an external environment. A MDP is defined as a 5-tuple {S, A, T, R, } , where S is a state space of the environment. A is an action space that can be taken by the agent to interact with the environment. T is a transition function defined as T(s, a, s � ) = p(s � |s, a) that tells a probability of next state s ′ if at state s, action a is taken. R is a reward function, i.e. R(s, a) , that returns a scalar value if the environment at state s and action a is taken. ∈ (0, 1] is a discount factor. The agent's task is to find an optimal deterministic policy ∶ S ↦ A , which is a mapping from states to actions, that maximizes an expected cumulative discounted return, denoted as J( ) a function of which is defined as where the expectation is with respect to the stochasticity of the transitions T . The term inside expectation is the infinite-horizon sum of discounted rewards. Alternatively, the policy can be stochastic and written as ∶ S × A ↦ [0, 1] that defines the probability of selecting actions at a state. Finding optimal policies MDP can be handled by the use of dynamic programming, e.g. value iteration or policy iteration algorithms Puterman [24].
Reinforcement learning is a family of learning algorithms in MDP in which J( ) is optimised w.r.t when the agent does not know the dynamics of the environment, i.e. unknown T and R . In unknown environment, planning techniques like value iteration or policy iteration can not be directly used. There are many different types of RL algorithms, but can often be categorized into two fashion: value-based and policy-based methods [32]. Value-based methods consist of Q-learning, SARSA, etc. that optimizes the policy by estimating intermediate terms, called value funcions or Q-value functions. Policy-based methods (direct policy search) optimize J( ) on the parametric policy space. The policy can be optimized using policy search techniques such as policy gradient [45], natural actor-critic [23], Bayesian policy gradient Ghavamzadeh et al. [7], Vien et al. [41], or black-box search such as CEM [20] and CMA-ES [15], etc. The first two techniques are local methods that can only look for locally optimal policies (more discussions over their differences can be found in [13]).
In the next sections, we first describe an application of CMA-ES for policy search as introduced in [15], then give a brief introduction of one policy gradient technique that models policies in reproducing kernel Hilbert space [19].

Direct policy search via CMA-ES
As direct policy search would directly search on the space of policies using stochastic optimization, one natural application of CMA-ES for direct policy search is to construct a search distribution on the policy space. Assuming that the policy space is parameterized by a function space H , in which each function h ∈ H is defined as h ∶ S ↦ ℜ . Specifically, each policy is written as a deterministic distribution as where ∈ A is an action, ∈ S is a state. For parametric policy approaches, the function space H is parametrized by a parameter space ∈ ℜ n . For example, h might be a neural network (where n is the number of trainable parameters in the network), or be a linear function of predefined features (where n is the number of features) as where j denotes the jth coordinate of , and j (s) is the jth feature. Therefore, the objective in Eq. 1 as a function of policy can be re-written as a function of parameters . The authors in [15] propose to use the CMA-ES algorithm in Algorithm 1 to optimize a cumulative reward objective function of , where the expectation is with respect to the stochasticity of the stochastic policy (s, a; ) and the dynamics T . For each candidate i , its evaluation in Step 5 in Algorithm 1 is J( i ) . Each J( i ) is computed using Monte-Carlo simulations. Specifically, a set of N trajectories are sampled by executing the policy ( i ) , then the evaluation is approximated as is a return of trajectory j, j . One advantage of CMA-ES is the ability to adapt the covariance matrix to control exploration, therefore it is considered to be more robust to noise and that is typically the case in robotics. The application of CMA-ES for direct policy search has shown many remarkable results in robotics Deisenrot et al. [4].

Policy gradient in reproducing kernel Hilbert space
An alternative direct policy search is policy gradient which uses gradient ascent to find a local optimal policy . In this section, we discuss how to use functional gradient ascent that would optimise the objective J(h) directly on a function space of h, instead of a parameterised function space through . We start by describing how to represent a function space in a non-parametric way through a reproducing kernel Hilbert space.
A reproducing kernel Hilbert space is a Hilbert space of functions in which the evaluation functional is a continuous linear functional. Denote S a set (e.g. a state space in MDP) and H a Hilbert space of functions on S . An evaluation functional is defined as a linear functional that has an evaluation at each point as Essentially, the above definition means the linear functional L evaluates f at a to receive L (f ) = f ( ).
According to the Riesz representation theorem [18], for all � ∈ S there exists a unique evaluation functional K( , ⋅) ∈ H satisfying the reproducing property: where ⟨⋅, ⋅⟩ denotes an inner product. The above result also means that we are evaluating the function K(⋅, � ) using the evaluation functional K( , ⋅) , and vice versa. We call the Hilbert space H with a reproducing kernel K, a reproducing kernel Hilbert space H K [29].
Recently, there have been efforts in integrating kernel methods with policy search in reinforcement learning. The work by [2,19] suggest to represent policies in reproducing kernel Hilbert spaces in order to enjoy rich representation and flexible complexity, in which each function h(⋅) in the definition of policies in Eq. 2 is a vector-valued function in RKHS H K , and defined as and K is a vector-valued kernel, K ∶ S × S ↦ L(A) [22], in which L(A) is the space of linear operators on A (an action space). In policy gradient, we often resort to stochastic policies in order to make gradients exist and promote exploration in RL, therefore Gaussian policies are commonly chosen and written as where is a covariance matrix that dictates the level of exploration. In all experiments, we use a diagonal matrix = 2 , where starts with a large value (when exploration dominates) and decreases gradually (until exploitation dominates) Lever and Stafford [19]. A more sophisticated strategy can be used to adapt , e.g. functional gradient Vien et al. [40].
Using this functional policy parametrization, we can analytically derive the functional gradient ∇ h J( h ) of the objective J( h ) w.r.t h by using the Fréchet derivative [18], and compute its approximate value given a set of sampled trajectories where R( ) = ∑ t t r t is the accumulated reward of the trajectory , and P( i ) is the probability distribution of trajectories. In the third step, we use the log trick is a sample of the Q-value function of ( t , t ) (this equality results from the Policy Gradient theorem [33]). Thus, the derivation in Eq. 7 results in ∇ h J( h ) as a function in H K as defined in Eq. 5. In addition, the functional gradient update at iteration k can be written as where is a step-size. As the representation of each policy h k+1 becomes more complex after each iteration [19] (the number of linear evaluation functionals K(s i , ⋅) involved in the representation of h k+1 ), a sparsification technique is often used to achieve a sparse and adaptive policy. A compatible kernel for Q-functions Q ( , ) can simply be derived based on kernel K, that is K h as Kernel regression methods can be used to approximate Q easily via kernel matching pursuit [42], kernelized Least Square Temporal Difference Q-learning (k-LSTDQ) [46], etc. This framework is called the compatible RKHS Actor-Critic framework (RKHS-AC) by [19]. Though RKHS-AC has very excellent policy modeling, it would only result in locally optimal policies. A global policy search technique like CMA-ES which may be considered as a powerful off-the-shelf tool but has many drawbacks on its current form. For many applications of CMA-ES, a parametric objective function is often designed, which requires a set of predefined features. Hence, the choice of features plays a very important role. Recently, the authors in [40] propose a natural actor-critic in RKHS algorithm by computing the Fisher information operator, which is analogue to the Fisher information matrix for parametric policies.

Problem statement
We consider a black-box functional optimisation problem [37] that finds the maximum of an unknown functional f ∶ H ↦ ℜ , where H = {h ∶ X ↦ Y} is a separable Hilbert space of vector-valued functions h with a domain X and output space Y . For each queried function h ∈ H , an evaluation y = f (h) is returned. The goal is to find an optimal function h * = arg max h∈H f (h) such that the functional f (h * ) is maximum with least queries, i.e. f (h * ) ≥ f (h), ∀h ∈ H . For example, we might want to find trajectories of minimum cost for a mobile robot to move from one location to a desired destination. A trajectory is represented by a function h ∶ [0, 1] ↦ C , that maps time t ∈ [0, 1] to robot configuration h(t) ∈ C where C denotes the configuration space of a robot. A configuration is the positions of all robot points as represented relatively in a fixed coordinate system. The cost functional f(h) maps each trajectory h to a scalar cost in ℜ . The cost might measure the efficiency of the trajectory (energy efficiency) and the proximity to obstacles (obstacle avoidance).

The CMA-ES-RKHS framework
We propose a new general-purpose CMA-ES-RKHS framework, a generalized CMA-ES variant in RKHS, that solves the above problem. We explicitly assume the function space H as a reproducing kernel Hilbert space (RKHS) associated with a kernel K. Each h ∈ H is defined as a mapping from an arbitrary space X to Y , h ∶ X ↦ Y . The function space H may be a vector-valued RKHS [22], denoted as H K , with the kernel K ∶ X × X ↦ L(Y) , where L(Y) is the space of linear operators on Y . For example, when X = ℜ n the simplest choice of K might be K(x, x � ) = (x, x � )I n , where I n is an n × n identity matrix, and is a scalar-valued kernel [29]. Each function h ∈ H is then represented approximately as a linear span of finite elements {x i , y i } as Now we define a search distribution over H . A direct extension of parametric CMA-ES is to use a Gaussian process over the solution function h, h ∼ GP(m, 2 C) , where m is a mean function in H K , C is a covariance operator on H K , and is a scalar global step-size. We discuss now how to update the functionals {m, C} and the parameter in our CMA-ES-RKHS framework, which is also summarised in Algorithm 2. Basically, one might think of the difference between CMA-ES versus CMA-ES-RKHS as: (1) the parameterization is a finite parameter space of ∈ ℜ n versus a potentially infinite space of h ∈ RKHS, (2) the search distribution is a multi-variate Normal distribution versus a Gaussian process, (3) each sample from the search distribution is a finite parameter vector i versus a function h i . We now describe the CMA-ES-RKHS framework in more detail.

Mean function update in RKHS
Assuming that at iteration k, we can sample a set of functions g i ∼ ℙ(0, C) (Step 4), where ℙ(0, C) is a Gaussian process with mean function 0 and covariance C. Many techniques for sampling from a Gaussian process are basically described in [25]. It is commonly known that a sample from a Gaussian process is not in H K with probability of 1, as discussed in detail by [1]. For any sampling techniques of a Gaussian process, we receive g i in a form of data tuples (x (i) , y (i) ) . In order to receive a function in H K , we approximate g i by a function g i ∈ H K . One technique allows us to do this is kernel ridge regression whose kernel is the covariance operator C(⋅, ⋅) (Step 5). Hence, in our framework each function g i is approximated by a function g i ∈ H K . As a result, a new function candidate sampled from the function distribution is h i = m + g i . The new mean function is updated as (Step 9), where the normalized weights w i satisfy As a result, after the update the new functional mean is an element in H K . We denote ḡ as There are a number of settings for w that might inherit from that of CMA-ES, for example: w i = 1∕ , w i ∝ − i + 1 ; or a better choice w i = log( + 1 2 ) − log(i) . In our experiment, we implement the latter as it performs consistently better in all domains.

Covariance operator update
The covariance operator update (Steps 14-15) is based on the best selected candidate functions in terms of their evaluations f (h i ) . Hence an empirical estimate of the covariance operator C on H K , called rank-update, is where c is a learning rate of rank-, and ⊗ denotes an outer product. Similar to parametric CMA-ES, we also consider the change of the mean function over time by estimating an evolution path function p c as (Step 14), where w is a variance-effectiveness constant, and c c is the backward time horizon for the evolution path function p c . This is low-pass filtered of chosen steps ḡ . As p c is just a linear combination of functions in H K , therefore p c is also an element in RKHS H K . As a result, a complete update of the covariance operator that combines both rank-1 and rank-is computed as (Step 15), where and c 1 , c are learning rates of rank-1 and rank-respectively. This reduces to a rank-1 update if c 1 = 1, c = 0 . Similarly, the update becomes a rank-update when c 1 = 0, c = 1.

Step-size update
The global step-size is adapted through the computation of a conjugate evolution path function p as (Step 11),ḡ where c is a backward time horizon for the conjugate evolution path function p . According to the bounded inverse theorem in functional analysis [3], C as computed in Eq. 13 is a linear operator in the RKHS H K , hence it has a bounded inverse C −1 . Therefore, p is updated in a way that renders it an element in H K . The volume and the correlation of the selected steps are compared to the expected value of the standard Gaussian process with a Dirac kernel. The fact that the former is larger than the latter makes increased, otherwise decreased. The update formula of (Step 12) is where ‖ ⋅ ‖ is the L 2 -norm, and d is a damping factor for . The term ‖GP(0, x )‖ is the expectation of all L 2 -norms of functions sampled from GP(0, x ) . This term can be computed in advance using Monte-Carlo simulations as where g i (⋅) is a function in H K approximated (via kernel ridge regression) from a sample g i drawn from GP(0, (⋅, ⋅)).

Sparsification and adaptive representation
We now discuss implementation concerns of the CMA-ES-RKHS algorithm. Firstly, the most critical one is the representation issue of mean functions m and covariance operators C. Secondly, it follows with discussions of parameter setting in CMA-ES-RKHS. Thirdly, we discuss how to deal with the update rule in Eq. 14 that involves finding the inverse operator C − 1 2 . Sparsification (Step 16): The updates of the mean function in Eq. 11 and covariance operator in Eq. 13 make their representation complexity increase linearly on the number of iterations. Though we receive an adaptive, flexible and complex policy, this would result in an expensive evaluation cost, e.g. computing m(x) or C(h, h � ) when needed, where x ∈ X and m, h, h � ∈ H K . Sparsification is a technique that is able to keep these kernel-based representation sparse and approximately accurate. Assuming that after each iteration, m and C are re-written in the forms of where N 1 , N 2 are the numbers of functions in the representation of m, C respectively, x i ∈ X, h i ∈ H K , and i , i ∈ Y . Assuming that a sparsification algorithm would sparsify m and C to become where x i ∈ X,h i ∈ H K and satisfying that n 1 ≪ N 1 , n 2 ≪ N 2 . In our CMA-ES-RKHS framework, we resort to two different sparsification techniques separately for m and C. We propose to use the kernel matching pursuit algorithm [42] to sparsify m. Its idea is to add kernel functions K(x i , ⋅) as features sequentially and greedily that maximally reduce the current approximation error. A tolerance constant is used to check the error reduction level before adding a new kernel feature. The use of tolerance also helps achieving more compact representation.
In general, we can use the kernel matching pursuit algorithm [42] to sparsify C. However, we aim to look for a method that will both sparsify C and together compute the inverse square root operator C − 1 2 , because C − 1 2 is also used in Step 11 in Algorithm 2. Therefore, we propose to use the kernel PCA method (kPCA) from [30] for achieving efficiently and fast both a sparse and compact covariance operator and its inverse square root operator. Specifically, we rewrite C in Eq. 16 as where H is a matrix whose ith column is h i , a N 2 × N 2 diagonal matrix = diag( i ) , and = H 1 2 . Via kPCA, C can be decomposed through a decomposition of the One could easily show that each eigenfunction of C is a linear span Hence v i is also an element in H K . From the decomposition of C via kPCA, we are now able to sparsify C by choosing only a small set of eigenfunctions of principal eigenvalues. Moreover, from the decomposition of C, the inverse square root operator C − 1 2 is derived as which is also a linear operator in RKHS H K . Parameter setting The mean function is represented by n 1 adaptive kernel features, hence it has n 1 pivotal parameters. In CMA-ES, parameter setting is based on the number of free parameters, the dimensionality of the search space. Though it is not precise, we use the same setting of CMA-ES for CMA-ES-RKHS's parameters, i.e. the parameters: c 1 , c c , c , c , d based on n 1 . We call n 1 the effective dimensionality on our CMA-ES-RKHS.

EDA and CEM in RKHS
For a short notice, the derivations of CMA-ES-RKHS can easily be transferred to derive EDA-RKHS and CEM-RKHS algorithms which are EDA and CEM in RKHS. Specifically, the updates at iteration k for EDA-RKHS and CEM-RKHS change step 9 and 15 in Algorithm 2 as follows Though the updates of C (k) EDA-RKHS and C (k) CEM-RKHS are simpler than C of CMA-ES-RKHS, they result in similar covariance operators on H K . Hence the implementation technique of EDA-RKHS and CEM-RKHS is similar to that of CMA-ES-RKHS as discussed above. Therefore, we will not put them in comparisons due to their weaker performance when compared to CMA-ES-RKHS.

Direct policy search via CMA-ES-RKHS
In RL literature, there is recent effort to model policies as functions in RKHS [2,19,40]. Such RKHS policy gradient approaches suffer from a problem of step-size and local optima. Though one of the extended work, RKHS EM-based policy search (RKHS-PoWER) by [40], would overcome this issue, it can only converge to local optima. We propose a new black-box direct policy search in RKHS that is based on CMA-ES-RKHS (the extension to EDA-RKHS and CEM-RKHS is similar). We use deterministic policies where each policy is a function in RKHS with a kernel K as first introduced in Eq. 5. Different from the standard CMA-ES based direct policy search [15], our non-parametric modeling enables optimisation in a very rich policy space and allows to learn more complex policies that is able to avoid local optima, enjoy adaptive and compact representation and do not depend on pre-defined features.
Adaptive CMA-ES direct policy search There is a naive way that modifies the parametric CMA-ES direct policy search [15] to become adaptive in selecting features. This naively proposed modification is used as a base-line to compare to CMA-ES-RKHS policy search. Assuming that we use a controller that is a linear span of RBF features i ( ) in which i is a center, where each i ( ) = K( , i ) Thus, h becomes a parametric function whose param- . This renders policies (a| ;w) a standard parametric policy. Therefore standard policy search algorithms like CMA-ES and policy gradient can be applied straightforwardly. Applying the CMA-ES based direct policy search method [15], the parameter space would be = {w i } ∈ ℜ n . We make a slight change to assume that the parameter space is {w i , i } , called adaptive CMA-ES direct policy search (CMA-ES-A). This adaptive algorithm would search for both optimum weights and optimum RBF features. A clear problem of CMA-ES-A is in determining the scaling of these parameters. However it is in principle overcome by the RKHS norm on functions which is similar to the setting in our CMA-ES-RKHS algorithm.

Experiments
In this section, we present experiment results to evaluate and compare our proposed framework to other state-of-the-art approaches. We first evaluate the advantages and general optimisation applications of CMA-ES-RKHS on two simple functional optimisation problems: 1-D and 2-D function spaces. We compare the behavior of CMA-ES-RKHS with other three other methods: the standard CMA-ES, the adaptive CMA-ES version (CMA-ES-A), and the functional gradient techniques. The next experiments are two RL tasks: inverted pendulum and double-pendulum. We compare our direct policy search via CMA-ES-RKHS to the standard CMA-ES policy search, the adaptive CMA-ES policy search, a parametric actor-critic, and the actor-critic in RKHS (RKHS-AC) methods. In all experiments, we use Gaussian Radial Basis Function (RBF) as an embedding kernel of the RKHS in which the bandwidths are set using median-trick. In other words, the bandwidth is set to the median (of all pair-wise distances between sampled data points) divided by the number of features used in each function. These experiments aim to evaluate the proposed CMA-ES-RKHS for: (1) the quality of the returned compact solution function, (2) the flexibility and power of the proposed method in capturing a complex solution function that can not be found easily by existing methods, (3) the applicability in practice, i.e. for direct policy search in reinforcement learning.

Synthetic domains
We design two unknown 1-dimension (1-D) and 2-dimension (2-D) functions f * . Each function is a mixture of two (multivariate in the case of 2-D) Gaussians, respectively. All optimizers are tasked to find a function h ∶ X ↦ ℜ , where h ∈ H K that minimizes the objective function as a square distance to the ground-truth. The objective is written as where x ∈ ℜ k , k = 1, 2 correspondingly to the 1-D or 2-D domain. We limit the domain X from x 0 to x T for 1-D task and in the box [x 0 , x T ] × [x 0 , x T ] for 2-D task. This task is a simplified version of many similar problems in machine learning and robotics, e.g. regularized risk functional (J(h) is a least-square cost and (x, f * (x)) are data samples) [29], trajectory optimisation (where f * (x) is a reference trajectory w.r.t time x ∈ [x 0 , x T ] ) [35], path planning or trajectory optimisation in RKHS [21], loss minimization inverse optimal control (where f * (x) is a demonstration trajectory) [5], etc. For example, in the case of path planning a robot must find a shortest path from one location to a destination while avoiding collisions with obstacles. The cost function in Eq. 19 is computing an Euclidean distance between the current configuration h(x) and the destination configuration represented as f * (x) , together with a distance to the closest obstacle. However, most of the above work must rely on discretization and parametric modelling. Functional gradient: Using functional gradient requires to know J and have access to the ground-truth function f * (CMA-ES-RKHS only accesses evaluations J(h)) from which we are able to use discretization to approximate the objective J in Eq. 19 as The functional gradient can be computed analytically as: Thus, a functional gradient update is h ← h + ∇ h J(h) . A sparsification technique [42] can be used to achieve a compact representation of h which renders the functional gradient approach an adaptive method too. That means the representation of h will be adapted to best approximate f * . Hence, discretization is required to be fine enough (T is large enough, we used T ≫ N ) to guarantee accurate approximation.

CMA-ES:
We assume that a parametric representation of h as a linear expansion of N features: We report the squared error J and the solution function for the 1-D task in Fig. 1. The results for the 2-D task are reported in Figs. 2 and 3. We create two versions for CMA-ES-A, one with good initialization (initial values of x t are centers for CMA-ES) and one with random initialization, called CMA-ES-A-R. In Fig. 1, the performance of CMA-ES-A-R is not good in terms of error. As demonstrated on the right picture, it can detect only one mode of the optimal function. Hence we stop reporting results from CMA-ES-A-R in other domains. One remarkable note is that CMA-ES initialization does not consist of two correct modes in its set of centers, hence it In the larger domain (2-D), CMA-ES-A performs much worse than our method. This is explained by the way our method approaches from a principled way, i.e kernel methods, for the scaling of parameters. The functional gradient method performs very well which also confirms that it can be very competitive when gradient information is known (in this case the form of J(h) is known). Fig. 2 shows very

Inverted pendulum
We use the same setting of the inverted pendulum domain as in [19]. This problem has an 1-dim action space that requires to bring the pole to the upright position and keep it balanced there. The dynamics of the system is � = + 0.05 + ; � = + 0.05a + 2 , is a small Gaussian noise N(0, 0.02 2 ) . We use N = 50 centers or features for all algorithms. We set = 0.99 and a horizon H = 400 . Each policy evaluation J(h) of policy h is averaged over 5 episodes.
The results of mean performance and its 95% confidence are computed over 15 runs and reported in Fig. 4. In this task, CMA-ES performs better than CMA-ES-A and CMA-ES-RKHS. CMA-ES-A has a much bigger search space comparing to that of CMA-ES, 3N versus N parameters. We conjecture that CMA-ES-RKHS performs worse because we use N as the effective dimensionality to set its parameters. Theoretically, CMA-ES-RKHS optimizes over a potentially infinite dimensional space. We tried to increase its effective dimensionality to 2N, called Eff. CMA-ES-RKHS. This modification improves the performance significantly. CMA-ES-A-R performs slowly but keeps improving constantly. Local direct policy search algorithms, AC and RKHS-AC, do not perform comparably to the other global direct methods. This experiment shows that CMA-ES-RKHS is able to avoid local optima.

Double pendulum
This problem consists of two links and two under-actuated joints. The system state is 4-dimensional of joint position and velocities = { 1 ,̇1, 2 ,̇2} . Actions are motor torques = [u 1 , u 2 ] , which are limited in [−5N, 5N] . The dynamics is simulated using second-order Runge-Kutta. We use a low sampling frequency of 50Hz at which torques could be applied. The start state is {0, − , 0, 0} . The reward function is r( ) = exp(−‖ − * ‖ W ) , where * = {0, 0, 0, 0} and W = diag(0.25, 0.0025, 0.25, 0.0025) . Each episode is simulated in 6s, which is equivalent to a horizon of 300 steps. Each policy evaluation J(h) of policy h is averaged over 5 episodes. We use N = 256 features or centers, and = 0.99 . The optimal policy returns 88. We only compare between global policy search methods via CMA-ES, CMA-ES-A, and CMA-ES-RKHS.
In this complex task, CMA-ES-RKHS has clearly outperformed other methods as seen in Fig. 5. Due to more expensive computation, we report an averaged performance over only three runs. CMA-ES performs worse because it still depends on fixed and pre-defined features, therefore in a more complex task it can only find a policy whose performance is up to the power and quality of the selected features. CMA-ES-A uses a non-principled scaling on its parameter space which inherently consists of two parts: the weights { } and the state information { i } . This scaling does not capture correctly distances between points on the parameter space, hence it leads to a non-optimal solution.

Conclusion
This paper proposes a CMA-ES-RKHS framework that generalises CMA-ES to handle functional optimisation in which the search is handled over a function space. The fact that the function space is modeled in reproducing kernel Hilbert space results in analytic update rules for CMA-ES-RKHS. On the other hand, the solution function attains compactness and flexibility characteristics. We apply CMA-ES-RKHS for direct policy search in which the policy is modeled in RKHS. Our experiments show that both CMA-ES-RKHS and direct policy search via CMA-ES-RKHS are able to represent a complex solution function compactly and adaptively. The result shows many interesting aspects and results of CMA-ES-RKHS: (1) explicitly handling functional optimisation in principle; (2) overcoming the issue of hand-designed feature in many practical applications of CMA-ES. Though offering many advantages, CMA-ES-RKHS is a kernel method therefore it also suffers from the problem of expensive computation. A study to investigate the way how to scale it will be a very promising research direction. There are also a number of other potential future research directions. A thorough study into hyperparameters of CMA-ES-RKHS will definitely be important to CMS-ES-RKHS. Moreover, more theoretical work and practical applications of CMA-ES-RKHS would be a very interesting future research direction.