Title: SHE 2 : STOCHASTIC HAMILTONIAN EXPLORATION AND EXPLOITATION FOR DERIVATIVE-FREE OPTI-MIZATION

Abstract: Derivative-free optimization (DFO) using trust region methods is frequently used for machine learning applications, such as (hyper-)parameter optimization without the derivatives of objective functions known. Inspired by the recent work in continuous-time minimizers, our work models the common trust region methods with the exploration-exploitation using a dynamical system coupling a pair of dynamical processes. While the first exploration process searches the minimum of the blackbox function through minimizing a time-evolving surrogation function, another exploitation process updates the surrogation function time-to-time using the points traversed by the exploration process. The efficiency of derivativefree optimization thus depends on how the two processes couple. In this paper, we propose a novel dynamical system, namely SHE 2 -Stochastic Hamiltonian Exploration and Exploitation, that surrogates the subregions of blackbox function using a time-evolving quadratic function, then explores and tracks the minimum of the quadratic functions using a fast-converging Hamiltonian system. The SHE 2 algorithm is later provided as a discrete-time numerical approximation to the system. To further accelerate optimization, we present P-SHE 2 that parallelizes multiple SHE 2 threads for concurrent exploration and exploitation. Experiment results based on a wide range of machine learning applications show that P-SHE 2 outperform a boarder range of derivative-free optimization algorithms with faster convergence speed under the same settings.

1. INTRODUCTION:

Derivative-free optimization (DFO) techniques (Powell, 1964) , such as Bayesian optimization algorithms (Snoek et al., 2012; Martinez-Cantin, 2014) , non-differentiable coordinate descent (Colson & Toint, 2002) , natural gradient method (Grosse & Salakhudinov, 2015; Hasenclever et al., 2017) , and natural evolution strategies (Schaul et al., 2011) , have been widely used for black-box function optimization. DFO techniques have been viewed as one of promising solutions, when the first-order/higher-order derivatives of the objective functions are not available. For example, to train large-scale machine learning models, parameter tuning is sometimes required. The problem to find the best parameters from the high-dimensional parameter space is frequently formalized as a black-box optimization problem, as the function that maps the specific parameter settings to the performance of models is not known (Golovin et al., 2017; Fazel et al., 2018; Xu et al., 2018b; Liu et al., 2018) . The evaluation of the black-box function is often computationally expensive, and there thus needs DFO algorithms to converge fast with global/local minimum guarantee.

Backgrounds. To ensure the performance of DFO algorithms, a series of pioneering work has been done (Conn et al., 2009; Powell, 2012; Jamieson et al., 2012; Augustin & Marzouk, 2017; Golovin et al., 2017) . Especially, Powell et al. (Powell, 1964; 1994) proposed Trust-Region methods that intends to "surrogate" the DFO solutions through exploring the minimum in the trust regions of the blackbox objective functions, where the trust regions are tightly approximated using model functions (e.g., quadratic functions or Gaussian process) via interpolation. Such two processes for exploration and exploitation are usually alternatively iterated, so as to pursue the global/local minimum (Augustin & Marzouk, 2017) . With exploration and exploitation (Dezza et al., 2017) , a wide range of algorithms have been proposed using trust region for DFO surrogation (Powell, 2002; Queipo et al., 2005; Snoek et al., 2012; Swersky et al., 2013; Schulman et al., 2015; Wu et al., 2017; Regier et al., 2017; Abdolmaleki et al., 2017; Osokin et al., 2017; Lyu et al., 2018; Arenz et al., 2018) .

Technical Challenges. Though trust region methods have been successfully used for derivative-free optimization for decades, the drawbacks of these methods are still significant:

• The computational and storage complexity for (convex) surrogates is extremely high. To approximate the trust regions of blackbox functions, quadratic functions (Powell, 2002; Regier et al., 2017) and Gaussian process (Snoek et al., 2012; Wu et al., 2017; Lyu et al., 2018) are frequently used as (convex) surrogates. However, fitting the quadratic functions and Gaussian process through interpolation is quite time-consuming with high sample complexity. For example, using quadratic functions as surrogates (i.e., approximation to the second-order Taylor's expansion) needs to estimate the gradient and inverse Hessian matrix (Powell, 2002; Regier et al., 2017) , where a large number of samples are required to avoid ill-conditioned inverse Hessian approximation; while the surrogate function in GP is nonconvex, which is even more sophisticated to optimize. • The convergence of trust region methods cannot be guaranteed for high-dimensional nonconvex DFO. Compared to the derivative-based algorithms such as stochastic gradient descent and accelerated gradient methods (Bottou, 2010; Su et al., 2014) , the convergence of DFO algorithms usually are not theoretically guaranteed. Jamieson et al. (Jamieson et al., 2012) provided the lower bound for algorithms based on boolean-based comparison of function evaluation. It shows that DFO algorithms can converge at Ω(1/ √ T ) rate in the best case (T refers to the total number of iterations), without assumptions on convexity and smoothness, even when the evaluation of black-box function is noisy.

Our Intuitions. To tackle the technical challenges, we are motivated to study novel trust region methods with following properties 1. Low-complexity Quadratic Surrogates with Limited Memory. To lower the computational complexity, we propose to use quadratic functions with identity Hessian matrices as surrogates. Rather than incorporating all evaluated samples in quadratic form approximation, our algorithm only works with the most-recently evaluated sample points. In this way, the memory consumption required can be further reduced. However, the use of identity Hessian matrices for quadratic form loses the information about the distribution (e.g., Fisher information or covariance (Hansen et al., 2003) ) of evaluated sample points. 2. Fast Quadratic Exploration with Stochastic Hamiltonian Dynamical Systems. Though it is difficult to improve the convergence rate of the DFO algorithms in general nonconvex settings with less oracle calls (i.e., times of function evaluation), one can make the exploration over the quadratic trust region even faster. Note that exploration requires to cover a trust region rather than running on the fastest path (e.g., the gradient flow (Hu & Li, 2017) ) towards the minimum of trust region. In this case, there needs an exploration mechanism traversing the whole quadratic trust region in a fast manner and (asymptotically) approaching to the minimum. Figure 1 illustrates the examples of exploration processes over the quadratic region via its gradient flows (i.e., gradient descent) or using Hamiltonian dynamics with gradients (Neal et al., 2011) as well as their stochastic variants with explicit perturbation, all in the same length of time. It shows that the stochastic Hamiltonian dynamics (shown in Figure 1(d) ) can well balance the needs of fast-approaching the minimum while sampling the quadratic region with its trajectories. Compared to the (stochastic) gradient flow, which leads to the convergence to the minimum in the fast manner, the stochastic Hamiltonian system are expected to well explore the quadratic trust region with the convergence kept.

Inspired by theoretical convergence consequences of Hamiltonian dynamics with Quadratic form (Su et al., 2014; Neal et al., 2011) , we propose to use stochastic Hamiltonian dynamical system for exploring the quadratic surrogates. 3. Multiple Quadratic Trust Regions with Parallel Exploration-Exploitation. Instead of using one quadratic cone as the surrogate, our method constructs the trust regions using multiple quadratic surrogates, where every surrogate is centered by one sample point. In this way, the information of multiple sample points can be still preserved. Further, to enjoy the speedup of parallel computation, the proposed method can be accelerated through exploring the minimum from multiple trust regions (using multiple Hamiltonian dynamical sys- Our work is inspired by the recent progress in the continuous-time convex minimizers (Su et al., 2014; Hu & Li, 2017; Xu et al., 2018a) on convex functions, where the optimization algorithms are considered as the discrete-time numerical approximation to some (stochastic) ordinary differential equations (ODEs) or dynamics, such as Itô processes for SGD algorithms (Hu & Li, 2017) or Hamiltonian systems for Nesterov's accelerated SGD (Su et al., 2014) . We intend to first study the new ODE and dynamical system as a continuous-time DFO minimizer that addresses above three research issues. With the new ODE, we aim at proposing the discrete-time approximation as the algorithms for black-box optimization.

3 2 1 0 1 2 X1 3 2 1 0 1 2 X2 (a) Gradient Flow 3 2 1 0 1 2 X1 3 2 1 0 1 2 X2 (b) Hamiltonian 3 2 1 0 1 2 X1 3 2 1 0 1 2 X2 (c) Stochastic GF 3 2 1 0 1 2 X1 3 2 1 0 1 2 X2 (d) Stochastic HAM

Our Contributions. Specifically, we make following contributions. (1) To address the three technical challenges, a continuous-time minimizer for derivative-free optimization based on a Hamiltonian system coupling two processes for exploration and exploitation respectively. (2) Based on the proposed dynamical system, an algorithm, namely SHE 2 -Stochastic Hamiltonian Exploration and Exploitation, as a discrete-time version of the proposed dynamical system, as well as P-SHE 2 that parallelizes SHE 2 for acceleration. (3) With the proposed algorithms, a series of experiments to evaluate SHE 2 and P-SHE 2 using real-world applications. The two algorithms outperform a wide range of DFO algorithms with better convergence. To the best of our knowledge, this work is the first to use a Hamiltonian system with coupled process for DFO algorithm design and analysis.

2. RELATED WORK AND PRELIMINARIES:

In this section, we first review the most relevant work of trust region methods for DFO problem, then present the preliminaries of this work.

2.1. TRUST REGION METHODS FOR DFO:

The trust region algorithms can be categorized by the model functions used for surrogates. Generally, there are two types of algorithms adopted: Gaussian Process (GP) (Snoek et al., 2012; Swersky et al., 2013; Wu et al., 2017; Lyu et al., 2018) or Quadratic functions (Powell, 2002; Regier et al., 2017; Osokin et al., 2017) for surrogation. Blessed by the power of Bayesian nonparameteric statistics, Gaussian process can well fit the trust regions, with confidence bounds measured, using samples evaluated by the blackbox function. However, the GP-based surrogation cannot work in high dimension and cannot scale-up with large number of samples. To solved this problem, GP-based surrogation algorithms using the kernel gradients (Wu et al., 2017) and mini-batch (Lyu et al., 2018) have been recently studied.

On the other hand, the quadratic surrogation (Powell, 2002) indeed approximates the trust region through interpolating the second-order Taylor expansion of the blackbox objective. With incoming points evaluated, there frequently needs to numerically estimate and adapt the inverse Hessian matrix and gradient vector, which is extremely time-consuming and sample-inefficiency (with sample complexity O(d 2 ) for d-dimensional DFO (Powell, 2002) ). Following such settings, (Regier et al., 2017) proposed to a second-order algorithm for blackbox variational inference based on quadratic surrogation, while (Osokin et al., 2017 ) leveraged a Gaussian Mixture Model (multiple quadratic surrogations) to fit the policy search space over blackbox probabilistic distribution for policy optimization. A novel convex model generalizing the quadratic surrogation has been recently proposed to characterize the loss for structured prediction (Arenz et al., 2018) .

In addition, some evolutionary strategies, such as Covariance Matrix Adaptation Evolution Strategy (CMA-ES) (Hansen et al., 2003; Abdolmaleki et al., 2017) , indeed behave as a sort of quadratic surrogate as well. Compared to the common quadratic surrogate, CMA-ES models the energy of blackbox function using a multivariate Gaussian distribution. For every iteration, CMA-ES draws a batch of multiple samples from the distribution, then statistically updates parameters of the distribution using the samples with blackbox evaluation. CMA-ES can be further accelerated with parallel blackbox function evaluation and has been used for hyperparameter optimization of deep learning (Loshchilov & Hutter, 2016) .

2.2. HAMILTONIAN SYSTEMS FOR CONVEX MINIMIZATION:

Here, we review the Nesterov's accelerated method for quadratic function minimization. We particularly are interested in the ODE of Nesterov's accelerated method and interpret behavior of the ODE as a Hamiltonian dynamical system. Corollary 1 (ODE of Nesterov's Accelerated Method). According to Su et al. (2014) , the discretetime numerical format of the Nesterov's accelerated method Nesterov (2004) can be viewed as an ODE as follow.

EQUATION

where f (X) is defined as the objective function for minimization and ∇ X f (Z(t)) refers to the gradient of the function on the point Z(t). Above ODE can converge with strongly theoretical consequences if the function f (X) is convex with some smoothness assumptions Su et al. (2014) .

Corollary 2 (Convergence of Eq 1 over Quadratic Loss). Let's set f (X) = 1 2 X -X * 2 2 According to the ODE analysis of Nestereov's accelerated method Su et al. (2014) , the ODE listed in Eq 1 converges with increasing time t at the following rate:

EQUATION

where X(0) refers to the initial status of the ODE. The proof has been given in Su et al. (2014) .

3. STOCHASTIC HAMILTONIAN SYSTEM FOR BLACK-BOX OPTIMIZATION:

In this section, we present the proposed Hamiltonian system for Black-Box minimization via exploration and exploitation. Then, we introduce the algorithms and analyze its approximation to the dynamical systems.

3.1. THE STOCHASTIC HAMILTONIAN EXPLORATION AND EXPLOITATION SYSTEM:

Given a black-box objective function f (X) and X ∈ R d , we propose to search the minimum of f (X) in the d-dimensional vector space R d , using a novel Hamiltonian system, derived from the ODE of Nesterov's accelerated method and Eq 1, yet without the derivative of f (X) needed. Definition 1 (Quadratic Loss Function). Given two d-dimensional vectors X and Y , we characterizes the Euclid distance between the two vectors using the function as follow.

EQUATION

)

where the partial derivative of Q on X should be ∂ ∂X Q(X, Y ) = X -Y indicating the fastest direction moving from X to Y .

Definition 2 (Stochastic Hamiltonian Exploration and Exploitation). As was shown in Eq. 4, a Hamiltonian system is designed with following two coupled processes: exploration process X(t) and exploitation process Y (t), where t refers to the searching time. These two processes are coupled withn each other. Specifically, the exploration process X(t) in Eq 4 uses a second order ODE to track the dynamic process Y (t), while the exploiting process Y (t) always memorizes the minimum point (i.e., X(τ )) that have been reached by X(t) from time 0 to t.

EQUATION

where (1) 

∂ ∂X Q(X(t), Y (t)) = X(t) -Y (t)

ζ t iid ∼ N (0, I) and ζ t ← ε • ζ t /|ζ t | 2 ; 6: X t ← X t-1 + α t • V t-1 ; 7: V t ← V t-1 + α t • Y t-1 + β t • X t-1 + γ t • X t + α t • ζ t ; 8: /* SHE 2 Y (t) Process Update*/ 9: if f (X t ) ≤ f (Y t-1 ) then 10: Y t ← X t ; 11: else 12: Y t ← Y t-1 ; 13:

end if 14: end for 15: return Y T ; ζ(t) 2 ≤ would help the system escape from an unstable stationary point in even shorter time. In the above dynamical system, we treat Y (t) as the minimizer of the black-box function f (X).

Indeed, SHE 2 approximates the black-box function f (X) using a simple yet effective quadratic function, then leverages the ODE listed in Eq 1 to approximate the minimum with the quadratic function. With the new trajectories traversed by Eq 1, the quadratic function would be updated. Through repeating such surrogation-approximation-updating procedures, the ODE continuously tracks the time-dependent evolution of quadratic loss functions and finally stops at a stationary point when the quadratic loss functions is no longer updated (even with new trajectories traversed). Remark 1. We can use the analytical results Su et al. (2014) to interpret the dynamical system (in Eq 4) as an adaptive perturbated dynamical system that intends to minimize the Euclid distance between X(t) and Y (t) at each time t. The memory complexity of this continuous-time minimizer is O(1), where a Markov process Y (t) is to used to memorize the status quo of local minimum during exploration and exploitation. Theorem 1 (Convergence of SHE 2 Dynamics). Let's denote x * as a possible local minimum point of the landscape function f (x). We have as t → ∞, with high probability, that X(t) → x * , where X(t) is the solution to (4).

Please refer to the Lemma 1 and Lemma 2 in the Appendix for the proof of above theorems. We will discuss the rate of convergence, when introducing SHE 2 algorithm as a discrete-time approximation to SHE 2 .

3.2 SHE 2 : ALGORITHM DESIGN AND ANALYSIS Given a black-box function f (x) and a sequence of non-negative step-size α t (t=0, 1, 2, . . . , T ), which is small enough, as well as the scale of perturbation ε, we propose to implement SHE 2 as Algorithm 1. The output of algorithm Y T refers to the value of Y t in the last iteration (i.e., the t th iteration). The whole algorithm only uses the evaluation of function f (x) for comparisons, without computing its derivatives. In each iteration, only the variable Y t is dedicated to memorize the local minimum in the sequence of X 1 , X 2 , . . . , X t . Thus the memory complexity of SHE 2 is O(1).

In terms of convergence, Jamieson et al Jamieson et al. (2012) provided an universal lower bound on the convergence rate of DFO based on the "boolean-valued" comparison of (noisy) function evaluation. SHE 2 should enjoy the same convergence rate Ω(1/ √ T ) without addressing any further assumptions.

3.3. APPROXIMATION ANALYSIS FOR SHE 2:

Here, we would demonstrate that the proposed algorithm behaves as a discrete-time approximation to the dynamical systems of X(t) and Y (t) addressed in Eq 4, while as α t → 0 the sequence of X t and Y t (for 1 ≤ t ≤ T ) would converge to the behavior of continuous-time minimizer -coupled processes X(t) and Y (t).

Given an appropriate constant step-size α t → 0 for t = 1, 2..., T , we can rewrite the the sequences X t described in lines 4-7 of Algorithm 1 as the following Stochastic Differential Equation (SDE) of X ω (t) with the random noise ω(t): 

EQUATION

EQUATION

The sequence Y t (t=0, 1, 2, • • • T ) always exploits the minimum point that has been already found by X t at time t. Thus, we can consider Y t is the discrete-time of Y (t) that exploits the minimum traversed by X(t). In this way, we can consider the coupled sequences of X t and Y t (for 1 ≤ t ≤ T ) as the discrete-time form of the proposed dynamical system with X(t) and Y (t).

4. P-SHE 2 : PARALLEL STOCHASTIC HAMILTONIAN EXPLORATION AND EXPLOITATION:

To enjoy the speedup of parallel computation, we propose a new Hamiltonian dynamical system with a set of ODEs that leverage multiple pairs of coupled processes for exploration and exploitation in parallel. Then, we present the algorithm design as a discrete-time approximation to the ODEs.

4.1 THE P-SHE 2 DYNAMICAL SYSTEM Given a black-box objective function f (X) and X ∈ R d , we propose to search the minimum of f (X) in the d-dimensional vector space R d , using following systems. Definition 3 (Parallel Stochastic Hamiltonian Exploration and Exploitation). As was shown in Eq. 7, a Hamiltonian system is designed with (1) N pairs of coupled exploration-exploitation processes: X i (t) and Y i (t) for 1 ≤ i ≤ N that explores and exploits the minimum in-parallel from N (random/unique) starting points, and (2) an overall exploitation process Y (t) memorizing the local minimum traversed by the all N pairs of coupled processes. Specifically, for each pair of coupled processes, a new surrogation model Q δ (X i (t), Y i (t), Y (t)) has been proposed to measure the joint distance from X i (t) to Y i (t) and Y (t) respectively, where δ > 0 refers to a trade-off factor weighted-averaging the two distances.

EQUATION

where

∂ ∂X Q δ (X i (t), Y i (t), Y (t)) = X i (t) -δ • Y i (t) -(1 -δ) • Y (t)

indicates the fastest direction to track Y i (t) and Y (t), jointly, from X(t). In the above dynamical system, we treat Y (t) as the minimizer of the black-box function f (X). Remark 2. We understand the dynamical system listed in Eq 7 as a perturbated dynamical system with multiple state variables, where all variables are coupled to search the minimum of f (X) through X i (t) (for 1 ≤ i ≤ N ). The memory complexity of this continuous-time minimizer is O(N ), where every Markov process Y i (t) is to used to memorize the status quo of local minimum traversed by the corresponding processes.

4.2 P-SHE 2 : SCALING-UP SHE 2 THROUGH PARALLELIZED FUNCTION EVALUATION The evaluation of black-box function f (x) is frequently time consuming. Algorithm 2 presents an algorithm namely P-SHE 2 that uses N SHE 2 threads to evaluate f (x) in parallel. Specifically, in each iteration, P-SHE 2 leverages the minimum among the N local minimums searched by the N parallel SHE 2 threads to accelerate the overall search. In P-SHE 2 , we use N threads of SHE 2 to search the minimum, where each thread uses three X j t , Y j t and V j t trajectories for exploration and exploitation. The memory complexity of P-SHE 2 is O(N ).

Algorithm 2 P-SHE 2 : Stochastic Hamiltonian Exploration and Exploitation 1: /* Starting P-SHE 2 by Initializing N instances of SHE 2 with Synchronization*/ 2: for j = 1, 2, 3, . . . , N do 3:

X j 0 , V j 0 ∼ randomly from R d and Y j 0 ← X j 0 ; 4: end for 5:

argmin ∀X∈{Y 1 0 ,Y 2 0 ,...,Y N 0 } f (X) → Y 0 ;

6: for t = 1, 2 . . . , T do 7:

for j = 1, 2, 3, . . . , N in Parallel do 8:

/* X(t) update for the j th SHE 2 thread*/ 9:

β t ← -(α t + 3/t) and γ t ← 3/t; 10:

ζ j t iid ∼ N (0, I) and ζ j t ← ε • ζ j t /|ζ j t | 2 ;

11:

X j t ← X j t-1 + α t V j t-1 ;

12:

W j t-1 ← δ • Y j t-1 + (1 -δ) • Y t-1 ;

13: As α t → 0, we can rewrite the sequences X j t , W j t and V j t described in lines 9-12 of Algorithm 2 as the following Ordinary Differential Equations (ODEs) of X j (t), W j (t) and V j (t):

V j t ← V j t-1 + α t • W j t-1 + β t • X j t-1 + γ t • X j t + α t •

EQUATION

where ζ j (t) refers to the continuous-time dynamics of sequence ζ j 1 , ζ j 2 , . . . , ζ j T and |ζ j (t)| 2 = ε for every time t. Through combining above three ODEs and Eq. 8, we can obtain the ODE of X(t) as:

EQUATION

Using same the settings, we can conclude that X j t would have similar behavior as X i (t) (for 1 ≤ i ≤ N in Eq 7). Thus, Algorithm 2 can be viewed as a discrete-time approximation of dynamical systems in Eq 7. Since the sequence Y t always exploits the minimum point that has been found by all N threads at every time t, we can use the algorithm output Y T as the minimizer of f (x). The proposed P-SHE 2 algorithm can be viewed as a particle swarm optimizer Kennedy (2011) with inverse-scale step-size settings. Compared to Particle Swarm, which usually adopts constant stepsize settings (i.e., α t , β t and γ t are fixed as a constant value), P-SHE 2 proposes to use a small α t , while setting β t = -(α t + 3/t) and γ t = 3/t for each (the t th ) iteration. Such settings help the optimizer approximates to the Nesterov's scheme, so as to enjoy faster convergence speed, under certain assumption. In terms of contribution, our research made as yet an rigorous analysis for Particle Swarm through linking it to to Nesterov's scheme Nesterov (2013) ; Su et al. (2014) .

5. EXPERIMENTS AND EMPIRICAL VALIDATION:

We provide three sets of experiments to validate our algorithms. In the first set of experiments, we demonstrate the performance of SHE 2 and P-SHE 2 to minimize two non-convex functions through the comparisons to a set of DFO optimizers, including Gaussian Process optimization algorithms (GP-UCB) (Martinez-Cantin, 2014), Powell's BOBYQA methods (Powell, 2009) , Limited Memory-BFGS-B (L-BFGS) (Zhu et al., 1997) , Covariance Matrix Adaptation Evolution Strategy (CMA-ES) (Hansen et al., 2003) , and Particle Swarm optimizer (PSO) Kennedy (2011) . For the second set of experiments, we use the same set of algorithms to train logistic regression (Lee et al., 2006) and support vector machine (Cortes & Vapnik, 1995) classifiers, on top of benchmark datasets, for supervised learning tasks. In the third set, we use P-SHE 2 to optimize the hyper-parameters of ResNet-50 for the performance tuning on Flower 102 and MIT Indoor 67 benchmark datasets under transfer learning settings.

5.1. NONCONVEX FUNCTION MINIMIZATION:

Figure 2 presents the performance comparison between P-SHE 2 , SHE 2 and the baseline algorithms using two 2D benchmark nonconvex functions-Franke's function and Peaks function. Figure 3 : Loss Minimization via DFO using Benchmark Datasets same computational/memory complexity, while P-SHE 2 (10) converges much faster than PSO(10). same phenomena can be also observed from the comparison between SHE 2 and PSO(1), both of which search with single thread. We can suggest that the adaptive step-size settings inherent from Nesterove's scheme accelerate the convergence speed. 

5.2. PARAMETER OPTIMIZATION FOR SUPERVISED LEARNING:

We use above algorithms to train logistic regression (LR) and SVM classifiers using Iris (4 features, 3 classes and 150 instances), Breast (32 features, 2 classes and 569 instances) and Wine (13 features, 3 classes and 178 instances) datasets. We treat the loss functions of logistic regression and SVM as black-box functions and parameters (e.g., projection vector β for logistic regression) as optimization outcomes. Note that the number of parameters for multi-class (#class ≥ 3) classification is #class × #f eatures, e.g., 39 for wine data. We don't include GP-UCB in the comparison, as it is extremely time-consuming to scale-up in high-dimensional settings.

Figure 3 demonstrates how loss function could be minimized by above algorithms with iterations by iterations. For both classifiers on all three datasets, P-SHE 2 (100)-the P-SHE 2 algorithms with N = 100 search threads, outperforms all rest algorithms with the most significant loss reduction and the best convergence performance. We also test the accuracy of trained classifiers using the testing datasets. Table . 1 shows that both classifiers trained by P-SHE 2 (100) enjoys the best accuracy among all above DFO algorithms and the accuracy is comparable to those trained using gradientbased optimizers. All above experiments are carried out under 10-folder cross-validation. Note that the accuracy of the classifiers trained by P-SHE 2 is closed to/or even better than some fine-tuned gradient-based solutions (Dogan & Tanrikulu, 2013; Roscher & Förstner, 2009) .

5.3. HYPERPARAMETER OPTIMIZATION FOR DEEP NEURAL NETWORKS:

To test the performance of P-SHE 2 for derivative-free optimization with noisy black-box function evaluation, We use P-SHE 2 to optimize the hyper-parameter of ResNet-50 networks for Flower (Nilsback & Zisserman, 2008) and MIT Indoor 67 classification (Quattoni & Torralba, 2009) tasks. The two networks are pre-trained using ImageNet (Krizhevsky et al., 2012) and Place365 datasets (Zhou et al., 2014) , respectively. Specifically, we design a black-box function to package the training procedure of the ResNet-50, where 12 continuous parameters, including the learning rate of procedure, type of optimizers (after simple discretization), the probabilistic distribution of image pre-processing operations for randomized data augmentation and so on, are interfaced as the input of the function while the validation loss of the network is returned as the output. We aim at searching the optimal parameters with the lowest validation loss. The experiments are all based on a Xeon E5 cluster with many available TitanX, M40x8, and 1080Ti GPUs.

Our experiments compare P-SHE 2 with a wide range of solvers and hyper-parameter tuning tools, including PSO, CMA-ES (Loshchilov & Hutter, 2016) , GP-UCB (Joy et al., 2016) and BOBYQA under the same pre-training/computing settings. Specifically, we adopt the vanilla implementation of GP-UCB and BOBYQA (with single search thread), while P-SHE 2 , PSO and CMA-ES are all with 10 search threads for parallel optimization. The experimental results show that all these algorithms can well optimize the hyer-parameters of ResNet for the better performance under the same settings, while P-SHE 2 has ever searched the hyperparameters with the lowest validation loss in our experiments (shown in Figure 4 ). Due to the restriction of PyBOBYQA API, we can only provide the function evaluation of the final solution obtained by BOBYQA as a flatline in Figure 4 . In fact, P-SHE 2 , PSO and CMA-ES may spend more GPU hours than GP-UCB and BOBYQA due to the parallel search. For the fair comparison, we also evaluate GP-UCB and BOBYQA with more than 100 iterations til the convergence, where GP-UCB can achieve 0.099854 validation error (which is comparable to the three parallel solvers) for Flower 102 task. Note that we only claim that P-SHE 2 can be used for hyerparameter optimization with decent performance. We don't intend to state that P-SHE 2 is the best for hyerparameter tuning, as the performance of the three parallel solvers are sometimes randon and indeed close to each other.

6. DISCUSSION AND CONCLUSION:

In this paper, we present SHE 2 and P-SHE 2 -two derivative-free optimization algorithms that leverage a Hamiltonian exploration and exploitation dynamical systems for black-box function optimization. Under mild condition SHE 2 algorithm behaves as a discrete-time approximation to a Nestereov's scheme ODE (Su et al., 2014) over the quadratic trust region of the blackbox function. Moreover, we propose P-SHE 2 to further accelerate the minimum search through parallelizing multiple SHE 2 -alike search threads with simple synchronization.

Compared to the existing trust region methods, P-SHE 2 uses multiple quadratic trust regions with multiple (coupled) stochastic Hamiltonian dynamics to accelerate the exploration-exploitation processes, while avoiding the needs of Hessian matrix estimation for quadratic function approximation. Instead of interpolating sampled points in one quadratic function, P-SHE 2 defacto constructs one quadratic surrogate (with identity Hessian) for each sampled point and leverages parallel search threads with parallel black-box function evaluation to boost the performance. Experiment results show that P-SHE 2 can compete a wide range of DFO algorithms to minimize nonconvex benchmark functions, train supervised learning models via parameter optimization, and fine-tune deep neural networks via hyperparameter optimization.

A PROOF OF THE CONVERGENCE OF SHE 2 DYNAMICAL SYSTEM Our goal is to show that in the system (4) we have X(t) → x * as t → ∞, where x * is a local minimum point of the landscape function f (x). Definition 4. We say the point x * is a local minimum point of the function f (x) if and only if f (x * ) ≤ f (x) for any x ∈ U (x * ), where U (x * ) which is any open neighborhood around the point x * .

Let us first remove the noise ζ(t) in our system (4). Thus we obtain the following deterministic dynamical system (X 0 (t), Y 0 (t)):

EQUATION

In the equations (a) and (b) of ( 10), the pair of processes (X 0 (t), Y 0 (t)) is a pair of coupled processes. In the next Lemma, we show that X 0 (t) converges to the minimum point of f (x) along the trajectory of X 0 (t) as t → ∞. Lemma 1. For the deterministic system (10), we have that X 0 (t) -Y 0 (t) 2 → 0 as t → ∞, so that the coupled process (X 0 (t), Y 0 (t)) converges to diagonal.

Proof. Set Ẋ0 (t) = P (t), we can write equation (a) in ( 10) in Hamiltonian form

EQUATION

Set H(X, Y, P ) = P 2 2 + Q(X, Y ). Then we have

EQUATION

As we have

EQUATION

)

Notice that by our construction part (b) of the coupled process (10), we have f (Y 0 (t)) ≤ f (X 0 (s)) for 0 ≤ s ≤ t. If X 0 (t) -Y 0 (t) = 0, then (X 0 (t) -Y 0 (t)) • Ẏ 0 (t) = 0. If X 0 (t) -Y 0 (t) = 0, then we see that Y 0 (t) = X 0 (s 0 ) for some 0 ≤ s 0 < t, and f (X 0 (s 0 )) < f (X 0 (t)). By continuity of the trajectory of X 0 (t) as well as the function f (x), we see that in this case Ẏ 0 (t) = 0, so that (X 0 (t) -Y 0 (t)) • Ẏ 0 (t) = 0. Thus we see that (13) actually gives d dt H(X 0 (t), Y 0 (t), P (t)) = -

EQUATION

)

From here, we know that H(X 0 (t), Y 0 (t), P (t)) keeps decaying until P (t) → 0 and X 0 (t) -Y 0 (t) 2 → 0, as desired.

Since f (Y 0 (t)) ≤ min 0≤s≤t f (X 0 (s)), Lemma 1 tells us that as t → ∞, the deterministic process X 0 (t) in (10) approaches the minimum of f along the trajectory traversed by itself. Let us now add the noise ζ(t) to part (a) of ( 10), so that we come back to our original system (4). We would like to argue that with the noise ζ(t), we actually have lim t→∞ min 0≤s≤t f (X(s)) = f (x * ), and thus X(t) → x * when t → ∞ as desired. Lemma 2. For the process X(t) in part (a) of the system (4), we have X(t) → x * as t → ∞, where x * is a local minimum of the landscape function f (x).

Proof. We first notice that we have ζ(t) = ε ξ t ξ t 2 , where ξ t ∼ N (0, I) is a sequence of i.i.d.

normal, so that ζ(t) 2 = ε, Eζ(t) = 0. Viewing (10) as a small random perturbation (see (Freidlin & Wentzell, 2012 , Chapter 2, Section 1)) of the system (4) we know that for any δ > 0 fixed, we have

EQUATION

as ε → 0. From here we know that the process X(t) behaves close to X 0 (t) with high probability, so that by Lemma 1 we know that with high probability we have lim Equation ( 19) enables the process X(t) to explore locally in an ergodic way its neighborhood points, so that if (17) is not valid, then X(t + dt) will explore a nearby point at which f (X(t + dt)) is less that min 0≤s≤t f (X(t)), and thus will move to that point. This leads to a further decay in the value of f (X(t)), which demonstrates that in he limit t → ∞ we must have (17), and the Lemma concludes.

Summarizing, we have the Theorem 1.

B DISCUSSIONS ABOUT THE RATE OF CONVERGENCE OF SHE 2 DYNAMICAL SYSTEM:

Here we provide a short discussion on the convergence rate of the algorithm SHE2. In the previous appendix we have demonstrated that the system (4) converges via two steps.

Step 1 in Lemma 1 shows that the differential equation modeling Nesterov's accelerated gradient descent (see Su et al. (2014) ) helps the process X(t) to "catch up" with the minimum point Y (t) on its path.

Step 2 in Lemma 2 shows that when t → ∞ the noise term ζ(t) helps the process X(t) to reach local

:



Figures and Tables:

Figure 1: Quadratic Surrogate Exploration via Gradient Flow, Hamiltonian Dynamics with Gradient, Stochastic Gradient Flow, and Stochastic Hamiltonian Dynamics with Gradient tems) concurrently while parallelizing the blackbox function evaluation to save the overall time. Certain synchronization among the Hamiltonian systems might be required to share knowledge among the parallel exploration-exploitation processes.

5) where ζ(t) refers to the continuous-time dynamics of sequence ζ 1 , ζ 2 , . . . , ζ T and |ζ(t)| 2 = ε for every time t. Through combining above two ODEs and Lemma 1, we can obtain the SDE of X(t) based on the perturbation ζ(t) as:

Figure 2: Performance Comparison using Nonconvex Functions 4.4 CONNECTION TO EXISTING SOLUTIONS

Figure  2.a and c present the landscape of these two functions, while Figure2.b and d present the performance evaluation of P-SHE 2 , SHE 2 and baseline algorithms on these two functions. All these algorithms are tuned with best parameters and evaluated for 20 times, while averaged performance is presented. Specifically, we illustrate how these algorithms would converge with increasing number of iterations. Obviously, on Franke's function, only P-SHE 2 (10), i.e., the P-SHE 2 algorithm with N = 10 search threads, CMA-ES, GP-UCB algorithms and BOBYQA converge to the global minimum, while the rest of algorithms, including SHE 2 and PSO, converge to the local minimum. Though P-SHE 2 (10) needs more iterations to converge to the global minimum, its per iteration time consumption is significantly lower than the other three convergeable algorithms (shown in Fig-ure2.e). The same comparison result can be also observed from the comparison using Peaks function (in Figures 2.b and 2.c, only CMA-ES and P-SHE 2 (10) converge to global minimum in the given number of iterations). Compared P-SHE 2 (10) to PSO(10), they both use 10 search threads with the

Figure 4: Hyper-Parameter Optimization for ResNet-50 on Benchmark Datasets 102 recognition(Nilsback & Zisserman, 2008) and MIT Indoor 67 classification(Quattoni & Torralba, 2009) tasks. The two networks are pre-trained using ImageNet(Krizhevsky et al., 2012) and Place365 datasets(Zhou et al., 2014), respectively. Specifically, we design a black-box function to package the training procedure of the ResNet-50, where 12 continuous parameters, including the learning rate of procedure, type of optimizers (after simple discretization), the probabilistic distribution of image pre-processing operations for randomized data augmentation and so on, are interfaced as the input of the function while the validation loss of the network is returned as the output. We aim at searching the optimal parameters with the lowest validation loss. The experiments are all based on a Xeon E5 cluster with many available TitanX, M40x8, and 1080Ti GPUs.

t→∞ f (X(t)) -min 0≤s≤t f (X(s)) = 0 .(16) Our next step is to improve the above asymptotic to X(t) → x * as t → ∞. Comparing with (16), we see that it suffices to showlim t→∞ min 0≤s≤t f (X(s)) = f (x * ) .(17)To demonstrate (17), we note that when t is large, we can ignore in (4) the damping term3 t), Y (t)) + ζ(t) = 0 , (b) : Y (t) = X(τ ) and τ = argmin τ ∈(0,t] , Y ) = X -Y we further see that the term ∂ ∂X Q(X, Y) also contribute little in (18). Thus part (a) of (18) reduces to a very simple equationẌ(t) + ζ(t) = 0 .(19)

Algorithm 1 SHE 2 : Stochastic Hamiltonian Exploration and Exploitation 1: X 0 , V 0 ∼ randomly from R d and Y 0 ← X 0 /* initialization by randomization*/; 2: for t = 1, 2, . . . , T do

Table content: <table><tr><td>3:</td><td>/* SHE 2 X(t) Process Update*/</td></tr><tr><td>4:</td><td>β t ← -(α t + 3/t) and γ t ← 3/t;</td></tr></table>

Accuracy Comparison for Parameter Optimization on Benchmark Datasets

Table content: <table><tr><td/><td>Iris</td><td/><td colspan="2">Wine</td><td colspan="2">Breast Cancer</td></tr><tr><td>Algorithms</td><td>LR</td><td>SVM</td><td>LR</td><td>SVM</td><td>LR</td><td>SVM</td></tr><tr><td>P-SHE 2 (100)</td><td>0.952 ± 0.023</td><td>0.987 ± 0.037</td><td>0.967 ± 0.034</td><td>0.961 ± 0.027</td><td>0.980 ± 0.014</td><td>0.982 ± 0.015</td></tr><tr><td/><td/><td colspan="3">Derivative-Free Solutions</td><td/><td/></tr><tr><td>SHE 2</td><td>0.493 ± 0.113</td><td>0.571 ± 0.065</td><td>0.628 ± 0.135</td><td>0.689 ± 0.138</td><td>0.848 ± 0.079</td><td>0.830 ± 0.112</td></tr><tr><td>PSO(100)</td><td>0.920 ± 0.059</td><td>0.933 ± 0.059</td><td>0.922 ± 0.047</td><td>0.900 ± 0.045</td><td>0.956 ± 0.028</td><td>0.947 ± 0.014</td></tr><tr><td>PSO(1)</td><td>0.361 ± 0.060</td><td>0.305 ± 0.034</td><td>0.355 ± 0.035</td><td>0.326 ± 0.015</td><td>0.608 ± 0.320</td><td>0.544 ± 0.211</td></tr><tr><td>BOBYQA</td><td>0.793 ± 0.118</td><td>0.933 ± 0.059</td><td>0.944 ± 0.042</td><td>0.953 ± 0.019</td><td>0.959 ± 0.022</td><td>0.869 ± 0.062</td></tr><tr><td>L-BFGS</td><td>0.747 ± 0.110</td><td>0.600 ± 0.202</td><td>0.772 ± 0.057</td><td>0.828 ± 0.093</td><td>0.815 ± 0.008</td><td>0.947 ± 0.014</td></tr><tr><td/><td/><td colspan="3">Derivative-based Solution</td><td/><td/></tr><tr><td>SGD*</td><td>0.950 ± 0.026</td><td>0.930 ± 0.038</td><td>0.975 ± 0.011</td><td>0.981 ± 0.016</td><td>0.965 ± 0.008</td><td>0.897 ± 0.041</td></tr></table>

Peer review (iclr2019):
 This paper suggests a continuous-time framework consisting of two coupled processes in order to perform derivative-free optimization. The first process optimizes a surrogate function, while the second process updates the surrogate function. This continuous-time process is then discretized in order to be run on various machine learning datasets. Overall, I think this is an interesting idea as competing methods do have high computational complexity costs. However, I’m not satisfied with the current state of the paper that does not properly discuss notions of complexity of their own method compared to existing methods. 1) “The computational and storage complexity for (convex) surrogates is extremely high.” The discussion in this paragraph is too superficial and not precise enough. a) First of all, the authors only discuss quadratic models but one can of course use linear models as well, see two references below (including work by Powell referenced there): Chapter 9 in Nocedal, J., & Wright, S. J. (2006). Numerical optimization 2nd. Conn, A. R., Scheinberg, K., & Vicente, L. N. (2009). Global convergence of general derivative-free trust-region algorithms to first-and second-order critical points. SIAM Journal on Optimization, 20(1), 387-415. I think this discussion should also be more precise, the authors claim the cost is extremely high but I would really expect a discussion comparing the complexity of this method with the complexity of their own approach. As discussed in Nocedal (reference above) the cost of each iteration with a linear model is O(n^3) instead of O(n^4) where n is the number of interpolation points. Perhaps this can also be improved with more recent developments, the authors should do a more thorough literature review. b) What is the complexity of the methods cited in the paper that rely on Gaussian processes? (including (Wu et al., 2017) and mini-batch (Lyu et al., 2018)). 2) “The convergence of trust region methods cannot be guaranteed for high-dimensional nonconvex DFO” Two remarks: a) This statement is incorrect as there are global convergence guarantees for derivative-free trust-region algorithms, see e.g. Conn, A. R., Scheinberg, K., & Vicente, L. N. (2009). Global convergence of general derivative-free trust-region algorithms to first-and second-order critical points. SIAM Journal on Optimization, 20(1), 387-415. In chapter 10, you will find global convergence guarantees for both first-order and second-order critical points. b) The authors seem to emphasize high-dimensional problems although the convergence guarantees above still apply. For high-order models, the dimension does have an effect, please elaborate on what specific comment you would like to make. Finally, can you comment on whether the lower bounds derived by Jamieson mentioned depend on the dimension. 3) Quadratic loss function The method developed by the authors rely on the use of a quadratic loss function. Can you comment on generalizing the results derived in the paper to more general loss functions? It seems that the computational complexity wouldn’t increase as much as existing DFO methods. Again, I think it would be interesting to give a more in-depth discussion of the complexity of your approach. 4) Convergence rate The authors used a perturbed variant of the second-order ODE defined in Su et al. 2014. The noise added to the ODE implies that the analysis derived in Su et al. 2014 does not apply as is. In order to deal with the noise the authors show that unbiased noise does not affect the asymptotic convergence. I think the authors could get strong non-asymptotic convergence results. In a nutshell, one could use tools from Ito calculus in order to bound the effect of the noise in the derivative of the Hamiltonian used in Lemma 1. See following references: Li, Q., Tai, C., et al. (2015). Stochastic modified equations and adaptive stochastic gradient algorithms. arXiv preprint arXiv:1511.06251. Krichene, W., Bayen, A., and Bartlett, P. L. (2015). Accelerated mirror descent in continuous and discrete time. In Advances in neural information processing systems, pages 2845–2853. Of course, the above works rely on the use of derivatives but as mentioned earlier, one should be able to rely on existing DFO results to prove convergence. If you check Chapter 2 in the book of Conn et al. (see reference above), you will see that linear interpolation schemes already offer some simple bounds on the distance between the true gradient of the gradient of the model (assuming Lipschitz continuity and differentiability). 5) Noise “The noise would help the system escape from an unstable stationary point in even shorter time” Please add a relevant citation. For isotropic noise, see Ge, R., Huang, F., Jin, C., and Yuan, Y. Escaping from saddle points-online stochastic gradient for tensor decomposition. Jin, C., Netrapalli, P., and Jordan, M. I. Accelerated gradient descent escapes saddle points faster than gradient descent. arXiv preprint arXiv:1711.10456, 6) Figure 2 Instead of having 2 separate plots for iteration numbers and time per iteration, why don’t you combine them to show the loss vs time. This would make it easier for the reader to see the combined effect. 7) Empirical evaluation a) There are not enough details provided to be able to reproduce the experiments. Reporting the range of the hyperparameters (Table 2 in the appendix) is not enough. How did you select the hyperparameters for each method? Especially step-size and batch-size which are critical for the performance of most algorithms. b) I have to admit that I am not extremely familiar with common experimental evaluations used for derivative-free methods but the datasets used in the paper seem to be rather small. Can you please justify the choice of these datasets, perhaps citing other recent papers that use similar datasets? 8) Connection to existing solutions The text is quite unclear but the authors seem to claim they establish a rigorous connection between their approach and particle swarm (“In terms of contribution, our research made as yet an rigorous analysis for Particle Swarm”). This however is not **rigorously** established and needs further explanation. The reference cited in the text (Kennedy 2011) does not appear to make any connection between particle swarm and accelerated gradient descent. Please elaborate. 9) SGD results Why are the results for SGD only reported in Table 1 and not in the figure? Some results for SGD are better than for P-SHE2 so why are you bolding the numbers for P-SHE2? It also seem surprising that SGD would achieve better results than the accelerated SGD method. What are the possible explanations? 10) Minor comments - Corollaries 1 and 2 should probably be named as theorems. They are not derived from any other theorem in the paper. They are also not Corollaries in Su et al. 2014. - Corollary 2 uses both X and Z. - Equation 5, the last equation with .dot{V}(t): there is a dot missing on top of the first X(t) “SHE2 should enjoy the same convergence rate Ω(1/T) without addressing any further assumptions” => What do you mean by “should”? - There are **many** typos in the text!! e.g. “the the”, “is to used”, “convergeable”,... please have someone else proofread your submission.