Efficient Hyper-parameter Optimization with Cubic Regularization
Authors: Zhenqian Shen
Hansi Yang
Yong Li
James Kwok
quanming yao
Venue: NeurIPS 2023 poster (2023)
Keywords: hyper-parameter optimization
cubic regularization

Abstract: As hyper-parameters are ubiquitous and can significantly affect the model performance, hyper-parameter optimization is extremely important in machine learning. In this paper, we consider a sub-class of hyper-parameter optimization problems, where the hyper-gradients are not available. Such problems frequently appear when the performance metric is non-differentiable or the hyper-parameter is not continuous. However, existing algorithms, like Bayesian optimization and reinforcement learning, often get trapped in local optimals with poor performance. To address the above limitations, we propose to use cubic regularization to accelerate convergence and avoid saddle points. First, we adopt stochastic relaxation, which allows obtaining gradient and Hessian information without hyper-gradients. Then, we exploit the rich curvature information by cubic regularization. Theoretically, we prove that the proposed method can converge to approximate second-order stationary points, and the convergence is also guaranteed when the lower-level problem is inexactly solved. Experiments on synthetic and real-world data demonstrate the effectiveness of our proposed method.

Efficient Hyper-parameter Optimization with Cubic Regularization
Zhengian Shen Hansi Yang²* Yong Li James Kwok² Quanming Yao
*Department of Electronic Engineering, Tsinghua University, Beijing, China Department of Computer Science and Engineering. Hong Kong University of Science and Technology, Hong Kong SAR, China {szg220mails, 1iyong070mail, qyaoaaλmail> tsinghua. edu. cn fhyangbw, jamesk) @cse. ust.hk
Abstract
As hyper-parameters are ubiquitous and can significantly affect the model perfor- mance, hyper-parameter optimization is extremely important in machine learning. In this paper, we consider a sub-class of hyper-parameter optimization problems, where the hyper-gradients the performance metric non-differentiable or the ous. However, existing algorithms, like Bayesian optimization and reinforcement learning, often get trapped in local optimals with poor performance. To address the above limitations, propose and avoid saddle points. First, we adopt stochastic relaxation, which allows obtain- ing gradient and Hessian information without hyper-gradients. Then, we exploit the rich curvature information by cubic regularization. Theoretically, we prove that the proposed method can converge approximate second-order stationary points, and the convergence is also guaranteed when the lower-level problem is inexactly solved. Experiments on synthetic and real-world data demonstrate the effectiveness of our proposed method.
not available. Such problems frequently appear when
is not continu-
cubic regularization to accelerate convergence
1 Introduction
In machine learning system, hyper-parameters are parameters that control the learning process. Well- known examples include learning rate, and choice loss functions or regularizers. Hyper-parametert play crucial role in determining the ultimate performance of machine learning models [1; 2: 3], affecting both the convergence of training algorithms and the generalization capabilities of the final models. Therefore, parrmmerrreeter optimization,
problem of finding the optimal values for the hyper-parameters in a machine learning system, is fundamental and important.
Since manual parameters, there have been many works on hyper-parameter optimization to automatically search for optimal hyper-parameters. Numerous existing srper--paramter optimization algorithms compute the hyper-gradient with respect hyper-parameter and utilize it to update current hyper-parameter [4; 6; 7]. However, computing the hyper- apper-gradient requires the objectives to be differentiable, and the pyperr-arameters to optimized must be continuous. As such, typer-gradients cannot be obtained in many important hyper- parameter optimization problems. For instance, in recommendation system and knowledge graph completion, non-differentiable evaluation metries such as MRR and NDCG are commonly used [8; 9; 10), which makes it impossible to compute hyper-gradients Moreover, many
tuning can be time-consuming and often inefficient to find better hyper-
"qqual contribution. Corresponding author.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

rrper-parameter' in machine learning procedures are discrete (e.g. number of layers and hidden dimension size), where computing their
is impossible either.
To solve these problems without the use of hyper-gradients, existing methods often use some derivative-free optimization algorithms, e.g., random search [3], Bayesian optimization [11] or genetic programming [12]. However, these methods are often very time-consuming, and cannot find optimal syper-parameters. Another type of method is based on stochastic relaxation [13; 14; 15], which constructs another differentiable objective from the original upper-level (i.e. the objective to find the aaper-parameter that have the best performance) objective, and use gradient- based optimization algorithms to optimize this new objective. This framework also covers methods based on reinforcement learning [16]. However, existing works under this framework only considers first-order information (i.e., gradient) of the new objective, which is still not sufficient enough as the upper-level objective can be very complex.
To address the limitations of existing works, in this paper, we propose a novel algorithm for hyper- parameter optimization that does not depend on hyper-gradients. Based on stochastic relaxation, we use cubic regularization [17] to optimize the relaxed objective and consider inexact solutions for the lower-level problems. Specially, our contributions can be summarized follows:
We propose to utilize cubic regularization based on stochastic relaxation to accelerate convergence and avoid saddle points pperrppareettrr optimization problems.
We provide theoretical analysis of the proposed method, showing that the proposed method converge approximate second-order stationary points with only inexactly solving lower level objectives.
Experiments on synthetic data and diverse applications demonstrate that the proposed method can help find hyper-pann saremmeresssrsss that lead models with good performances.
2
Related Works
2.1 prper-parameter Optimization
Plyper-parameee optimization [18; 19], as a subfield of automated machine learning (AutoML) [20]. aims automate the design and tuning of machine learning models. Most existing works on hyper- parameter optimization adopt the bi-level formulation [21; 4], where hyper-parameters are optimized the upper-level objective, and model parameters are trained in the lower-level objective.
Generally, hyper-paramett optimization algorithms can be divided into two types. The first is derivative-free methods, which only depends the objective value. For example, grid search [22] and random search [3] are two simple search methods to try different choose the best among them based the upper-level objective. Another example is model-based Bayesian optimization [11; 2], which builds a probabilistic model to map different hyper-parameters their evaluation performance and updates this model along with the optimization process. Genetic programming [12] starts with population of randomly sampled
values and
and evolves new generations of syper-parameters by selecting top-performing ones and recombining them. Hyperband [23] improves upon random search by stopping training models with bad performances early before convergence. These methods may suffer from huge computation cost when the number of hyper-parameters is more than a few, and usually do not have theoretical guarantee for convergence.
Another approach uses hyper-gradient update hyper-parameters by gradient descent [4; 5; 6; 7]. Despite their solid theoretical analysis and good empirical performances, these methods require both continuous hyper-parameter and differentiable upper-level objective compute the hyper-gradient. In contrast, our method handles the cases where hyper-gradient is not available, such as when the parametrssssrss are discrete or the objectives are dondifferrntiable. Moreover, exact computation of the grppiet can be cine-consuming, and is only feasible for small models.
2.2 Cubic Regularization
Cubie regularization [17] a novel technique for general unconstrained and smooth optimization problems. The main idea is to construct cubic polynomial from the original problem, and update the parameters by finding the minimizer of this polynomial in each iteration. Since the pioneering work [17] in this direction, many variants have been proposed to improve its performance. For

---

example, ARC [24] proposes to use approximate global minimizer to reduce the computation cost of cubic regularization method. More recently, cubic regularization is utilized to accelerate convergence in stochastic non-convex optimization problems [25; 26; 27].
Compared to first-order optimization methods, cubic regularization escape saddle points [26]. Despite its solid theoretical foundations and good empirical performances, cubic regularization involves computing the full Hessian matrix, which prevents its application when the number of parameters is large (e.g., training a neural network). Nevertheless, since hyper- parameters are often -iwdimensional, computing the Hessian matrix will not be bottleneck. Our proposed method is also the first to use cubic regularization in a bi-level optimization problem.
converge faster [17; 25] and
3 Methodology
3.1 Problem Formulation
optimization with stochastic relaxation [13] can be formulated as the following bi-level optimization problem:
Mathematically,
0- = arg min [(()) (((0)=Bsmpoff Eg-yo (=) [ (x),, (3) Prail] } s.t. w" w²(=) (=) arg min F(w, Dua); (1)
0
where w denotes the model parameter and z denotes the hyper-parameter. Through stochastic
relaxation, the hyper-parameter 2 is sampled from the probability distribution po(z) parameterized
by 0. H and F are the performance evaluation metrics for the training dataset Dum and the validation dataset Deal, respectively.
Different from directly searching for good hyppe--parmmters, (I) introduces probability distribution (2(2) on hyper-paramede parameter 2, so that we can optimize 0 instead of 2. The following Lemma shows the validation set:
that optimizing 0 is equivalent to finding the best hyper- paremmeer
Lemma 3.1. Assume that for any z", there exists O (could be infinite) such that po(z) ō(= x²) where 8() is the Dirac's 's delta function. Then, (i) inf I (0) = min, H(w" (z), Dvai); (ii) The optimal 0 in (1) satisfies pz(2) 8(2 arg min, H(w (2), 2; Dval))-
The proof is in Appendix A. The solution of (1) can be easier as [(0) in (1) is always differentiable, even when we cannot compute the hyper-gradient with respect 2. As introduced in Section 1, hyper-parameter optimization problems where hyper-gradient cannot be computed are really common. Here are two examples:
Example 1: Hyper-parametio optimization for knowledge graph completion [9; 28J. A knowledge graph G is set of triplets {(h, r, t)} that indicates entity h has relation another entity t. Knowledge graph completion is to learn a model to predict unknown facts (i.e., triplets with missing components). Here, we mainly focus embedding-based models, and consider the task of tuning hyper-parameters for this type of model. Under the formulation of (1), w refers to the embedding model parameters, refers the hyper-parametent (e.g. learning rate, regularizer, scoring function). И is discrete evaluation metric (c.g. MRR, @it@ 10). F' corresponds to the loss function for training. Dua and Dval are formed by training and validation triplets in the knowledge graph.
Example 2: Schedule search for learning with noisy labels [29; 30J. In this application, we have a noisy training dataset (i.e., labels may be incorrect), and we use proportion schedule control the proportion of small-loss samples for model training in each epoch, which can avoid over-fitting on incorrect labels. Existing works [29; 30] demonstrate that different schedules have significant impact on the final performance. Therefore, we consider using hyper-paramen and we can optimize
to express this schedule, help find good schedules that lead good models. Under the formulation of (1), we let to be the model parameter, w" be the final model parameters trained with schedule parameterized by 2. F denotes the training loss and H denotes the negative validation set accuracy. Dma is the training dataset with noisy label, and Deal is the clean validation dataset.
3.2 Proposed Method
Since the upper-level objective in (1) can be highly non-convex, common optimization methods such as gradient descent and Newton method may get stuck in saddle points, which leads to poor
3

---

performance. As such, we propose to use cubic regularization, which can avoid convergence to saddle points, to optimize the relaxed objective J I()) in (1). We first briefly introduce the main procedure of using cubic regularization to minimize given objective J [0).. Denote tth as the value of & in the m-th iteration, cubic regularization will first construct the following surrogate:
Cy A T(mm)) (VJ()'))
(VTT()") A P
3
(2)
where is a parameter set in advance, We will then update n by 0 n + arg mina Cm (A) which involves minimizing the above surrogate as sub-procedure.
To compute gradient V (O")) and Hessian (OT),, we have the following propositions:
Proposition 3.2. V (0) = pe (8* (0; 20). where 8" (0; 2) = [1ww Dwv)) V VoogPo(z).
Proposition 3.3. (2) = Bpo(=) (0; where H*;0; = [z)" Dvat)
(V (o)) log
2 po
Proofs of these two propositions are in Appendix A. Note that do not require the original upper- level objective H to be differentiable, as we only need to use the value of upper-level objective H (ww (2), Dvat) to compute the gradient and Hessian of [(). Nevertheless, these two propositions still cannot be directly used, as they require infinite number of samples z as well as the corresponding exact lower-level solutions (z). and both of them are impossible obtain in practice. As such, we first define (2) (2) to be an approximate solution of the lower-level objective with hyper-parameter H(w' (z), 2; D v lo (z)) approximations of 8" (0; and (0; 2). With these
2. Then we define s0; 2) (ww (z), 2; Dval) Vlog and H0:
(V log (=) V log pel notations, the approximated gradient (denoted as gnt) and Hessian (denoted as 3") in the m-th iteration are computed as:
K E
K Σ K
1
g- -"=
33)
B"D
(4)
k=1
In other words, we will draw K different hyper-parameters 2 from the distribution pam (2) in m-th iteration, and use them to obtain the corresponding approximate solutions w w'(*)). Then we can compute the approximate gradient g" and Hessian B" and approximate surrogate En (0) is then given by:
2m = (g A A D
A
(5)
Compared to (0) in (2), we remove the constant term J(0"), which will not affect the minimizer.
is then updated as gm-1 9m + m where m. arg min غ. (A) is the minimizer of this new surrogate and can be obtained by existing solvers on cubic regularization (e.g., [24]). The
complete algorithm shown in Algorithm
Algorithm 1
optimization with cubic regularization.
1: Initialize 0 and =* randomly sampled from pgo (z).
2: for 0,, ... M - 1 do
3: for == 0.... K do
4: Draw hyper-parameter R from Pa (2);
5: Optimize lower-level objective in (1) to obtain w' v=(=))) most expensive step
6:
Update 2n = ² if we have better validation performance with hyper-parameter &
7:
end for
8:
Compute g" using (3) and 8m using (4):
9: Compute A m by m arg mign &m (A):
10: Update 0m+ = mm m:
end for
12: Perform the training step with x¹ and obtain the final model parameter w":
13: return Final model parameter w*
Since the number of hyper-parameters is often not too large, also low-dimensional variable in most cases. For example, for the knowledge graph application in Section 4.2. the dimension of

---

0 is only 25, as we sum up the number of possible values for each hyper pareer-rrrmtter. This is almost negligible compared with the embedding model, which often has millions of parameters. As such, optimizing (5) takes little time compared with solving the lower-level problem, which will also be empirically verified in Section 4.4.
3.3 Theoretical Analysis
Previous works that use cubic regularization for stochastic non-convex optimization [23; 26] do not consider bi-level optimization problems, and assume unbiased estimations of gradient and Hessian. As such, their analysis cannot be directly generalized here approximated gradient and Hessian are not unbiased due to inexact lower-level solutions. To analyze the convergence Algorithm 1, we first introduce the definition for saddle point and approximate second-order stationary point, which is also used in [26] for convergence analysis:
Definition 3.4 ([26)). For an objective I with p-Lipschitz Hessian, 0 is an e-second-order statioon point of J if J l and Amin (V²J(O) -ype.
The definition of c-second-order stationary point can be regarded as a extension of c-first-order stationary point. For e-first-order stationary points, we only require the gradient (O) to satisfy [V)]()) S e. And for e-second-order stationary points, we also require the eigenvalues of its Hessian be large enough, so as avoid getting stuck in saddle points in Definition 3.5.
Definition 3.5. For an objective I with p-Lipschitz Hessian, is e-saddle point of J if (O)]()) S e and Amin ( (O) < <vp..
Now we make the following assumptions. Assumption 3.6 (i) can be easily satisfied as the upper- level objective H is naturally should match (at least not smaller than) the upper bound p in Assumption 3.6 (i). This can be done
in real-world applications, and the parameter p in (2)
by gradually increasing p in (2) to ensure convergence. Assumption 3.6 (ii) is also satisfied by common probability distributions used in our experiments (e.g- sigmoid function for experiments on synthetic data and softmax liked distribution for hyper-pane parameter optimization for knowledge graph completion). Assumption 3.6 (i), (iii) and (iv) are popularly used in existing works on cubic regularization [25; 26).
Assumption 3.6. Following assumptions on I are made:
(i) Objective: infg (0) > oo and J is second-order differentiable with p-Lipschitz Hessian.
iii Bounded gradient and Hessian probability distribution function (= V log (o (2) s Q and ||V log (e(2)) 2|2 < Q2.
(iii) Bounded variance and bias in gradient estimation: (x) (0) 2) (8 (s(; of and (0s =) (s) 2)) < M, for all 0 and z.
(iv) Bounded error Hessian estimation: denote the spectral norm of matrix (i.e., the maximum singular value of A) as ||sp we have Eps(z) [[; - Epo(s) H ( - E (x) (0; lsp S M2 for all and 2.
< 02 and
The following theorem guarantees the convergence of Algorithm (proof is in Appendix B).
Theorem 3.7. For any J satisfying Assumption 3.6, and approximate lower-level solution w' w'(z) with (w" (z), 2; Dwa) H (w' (w'(2), Dwat)l S min ( (हáq,) +2();+Q2)> there exists K = ( log(t)), with probability 1 for any 8 > 0, such that Algorithm produces an e-second-order stationary point M of J, with M = 0 (0)) -3((0))
w/p
2.5
Theorem 3.7 shows that cubic regularization (Algorithm 1) takes only M ((1/e1.5) iterations and KM = (1)c2-5) samples to obtain an e-second-order stationary point. This is strictly faster than gradient-based optimization methods [14; 15), which require (11/e samples only obtain an
--first-order stationary point. We also note that the convergence rate depends on two parts: inexactness
'In other words, for any z, n, have [=) V>J(i))1 I e ||= vllg- where Il denotes the matrix spectral norm.
5

---

of lower-level solutions [(w 2; Dvai) - (2) Dra)). and the number of hyper parameter samples K to compute the approximate gradient g"' and Hessian BT.. To better demonstrate their impacts, we consider the following two special cases:
Set K oe: in such case, the convergence is solely controlled by the inexactness of lower-level solutions. To obtain a good & (hence good syper-parameter 2), we need to obtain good lower- level solutions, which intuitive as inaccurate lower-level solutions cannot help us identify good hyper-parameters
Using exact lower-level solutions: in such case, the convergence solely controlled by the number of syper-paramete samples K. To obtain good 0 (hence good be sufficiently large to accurately compute the approximate gradient g" and Hessian BM. This is also intuitive as inaccurate gradient or Hessian cannot lead to good 0, and matches previous analysis [25; 26] that focus on single-level optimization problems.
2), we need K to
The above analysis will also be justified by empirical results in Section 4.3, where we consider the impact of different number of samples K and different training epochs in lower-level, which causes the inexact solution of the lower-level objective.
Comparison with existing works While there have been many works on hyper-parameter opti- mization [30; 14; 15] based on stochastic relaxation [13], our method has two significant differences. First, the best of our knowledge, our method is the first method with convergence guarantee to approximate second-orde, stationary points (second-order convergence guarantees), while previous methods only have first-order convergence guarantees [30; 15]. As such, our method can avoid convergence to bad saddle points. While recent work [31] also claims that they are able to escape saddle points in bi-level optimization, they require
to be available and more strict assumptions on objectives. Moreover, our method allows the use of inexact lower-level solutions, while existing methods all assume lower-level solutions are exact. Since lower-level problems are often complex non-convex problems, inexact lower-level solutions are more realistic and suitable for convergence analysis. Detailed comparison upon existing works can also be found Appendix C.
Experiments
4.1 Experiment Synthetic Data
First, we consider problem of filtering useful features for linear model [32; 33]. We construct a synthetic dataset on 5-way classification, where the inputs are 50-dimensional vectors. Of the 50 features, 25 have different distributions based on their ground-truti labels, generated from Gaussian distributions with different means for each class and the same variance. The remaining 25 features are filled with i.i.d. Gaussian white noise. The vector to specify which dimensions are masked. We use cross entropy loss as the loss function in the lower-level objective and use AUC to measure the performance of the classifier in the upper-level objective. For the mask of i-th dimension 2;, we use sigmoid function to represent the probability to mask that dimension, i.c., Pe, (z: = 1) = 1/(1+exp(-9:)) and Pe, (z, = 0) = 1 po. (z+ 1). The complete distribution (=(=) is represented by successive multiplications pa
& {0. 1) )50 is 50-dimensional 0/1
1 Po, (z:).
Experiments conducted on a 24GB NVIDIA GeForce RTX 3090 GPU.
We compare the proposed method (Cubic) with methods that are also based on stochastic relaxation, including gradient descent (GD) [34], natural gradient (NG) [35], and Newton's methods (New- ton) [36]. The validation AUC of different methods are shown in Figure ((a). We can see that our proposed method clearly outperforms other methods based stochastic relaxation.
Furthermore, to empirically verify that our proposed method does escape bad saddle points, we compare the eigenvalues of Hessian matrix (²) with 0 found by different methods. For easier comparison, we divide each method's eigenvalues by their largest, which makes different methods have the same largest eigenvalue as 1. As can be seen from Figure ((b), only our proposed method obtains an optimized 0 that make the negative eigenvalues less significant compared to positive eigenvalues, which means only our proposed method can better escape the saddle points in J (()).

---

1.00
Q
036
MO
0.75
Nanator
0.92
velteattoon AUK.
ReMtiv egennttis
0.25
0.00
G0
BFG5 Outic
Neaton
-0.75
100
20c
400
500
000
10
20
30
50
Running timett
Number of elgervaizes
(a) Validation AUC comparison.
(b) Eigenvalue comparison.
Figure 1: Comparison among different gradient based methods for mask learning in synthetic data
classification.
4.2 Experiments on Real-world Data
In this section, we mainly consider two applications, hyper- hyper-parameter optimization for knowledge graph completion and schedule search on learning with noisy training data, which are both introduced as examples Section 3.1. For baselines, we choose hyper-parameter optimization methods that are commonly used in hyper-parameter optimization literature [19; 30; 28), including random search (Random) [3], Bayesian optimization (BO) [11], Hyperband [23] and reinforcement learning (RL) [37]. These baseline methods applied to the original hyper-pate pprrmeee search problem instead of the problem after stochastic relaxation.
4.2.1 Hyper-parameter Optimization for Knowledge Graph Completion
For this application, need tune several discrete hyper-paramene for knowledge graph embedding model. The hyper-parameters to be tuned include negative sampling number, regularizer, loss function, gamma (for margin ranking loss), initializer, learning and score function, and more details can be found in Appendix E.2.
As all ayper-parameters. considered in this application discrete, choose softmax-like -iked distribu- tions to represent the probability to select specific value for each aaper-parameter For example, consider haperr-praaetter 2; with possible values {ai1...., J. the probability of 2ij is
p(2; = -iy) Σe= "xp( where O,j denote the probability function parameter corresponding to hyper-parameter Zij. The dimension of distribution parameter is 25. Two well-known knowledge graph datasets, FB 5k2337 [38] and WN18RR [39], are used in experiments and their statistics are in Appendix E.2.
exp(O.))
The experimental results are presented in Figure 2, which shows the achieved maximum reciprocal rank (MRR)* 4 with respect to the number of trained models for different search algorithms. Our proposed method exhibits rapid convergence, as Theorem 3.7 guarantees its fast convergence rate. Moreover, our proposed method outperforms other search algorithms and is more close to the global optimum in terms of MRR because cubic regularization allows our proposed method to escape from bad saddle points, while other methods may get trapped in such saddle points.
4.2.2 Schedule Search on Learning with Noisy Training Data
For this application, we consider the problem of searching for schedule (t) parameterized by
hyper-parameter [30]. The typer-parameter is divided into two parts: 2 = (a. { }), and R2 (t) parameterized as follows: (t El= arr (t; B.), where a (aj..... a,) can be seen as set of weights with a; 2 0 and Σ a; 1. Settings of functions r1: 3 parameterized with 3
shown Table 1. Here the dimension of hyper-parameter 36.
Denote the number of test triplets as n, and the rank of mtth triplet prediction as rank... MRR

---

0.32
Eest Nestvalettion
0.400
0.335
Bandon
0.350
Random
0.24
80
50
Iyperband
0.325
Hyporband
6.22
Cubic
0.300
Cutee
2000
4000p
54000
10000
20000.
50000
42000
50000
60000
70000
80000
Burning timess
Running Orners
(a) FB15K237.
WN18RR.
Figure 2: Validation MRR w.r.t. total number of trained models during search in optimization for knowledge graph completion.
We use CIFAR- 10 dataset for experiments with 50k, 5k, 5k image samples as training, validation, test data, respectively. We consider two different label noise settings: 50% symmetric noise and 45% pair flipping noise, and detailed generation process for the two settings is introduced in Appendix E.3.
Figure 3 shows the validation accuracy with respect the total number of trained models in the search process. As can be seen, our method with cubic regularization is more efficient than the other algorithms compared. Moreover, our proposed method converges a more accurate model because of the ability of cubic regularization to escape from bad saddle points.
90
boot
Bendem
Best
37.5
Random
35
35.0
Hyper Cubic
Hyoe fuand
32.5
Cuolc
20000 00000
atdno 100000 120000
Running 140000
Running
(a) symmetry flipping (50%).
(b) pair flipping (45%).
Figure 3: Validation accuracy w.r.t. total number of trained models during search for schedule search on learning with noisy training data.
4.3 Ablation Study
Table 1: The basis functions used to define the search space in the experiments. Here, I' is the total number of training epochs.
Here conduct ablation study two important fac- tors of our proposed method: the number of samples K (which affects the error in estimating the gradient and Hessian of J(0)) and number of training epochs for each model (which causes the inexact solution of the lower-level objective). The experiments are done on the synthetic data experiment section 4.1. Results are shown in Figure 4, where in the first row we fix the number of upper-level iterations among different parameter settings, while overall training time is con- trolled the same the second row. We fix the epoch number as 40 for experiment on K and set K = 5 for experiment on epoch number. In the experiment on the number of samples K (Figure 4(a)), we can see that our method performs better when K is larger, which means that when the estimation for gradient and Hessian of J a()) the more precise, the optimization of J(0) is more effective. However, when the overall training time is limited, excessively large K (c.g. K 50) causes decline
(1:
38 (+)
93188
333 log( 10
(1;
(1;
398 (7() 00
(1+236))
102
(tt
34
(1+821))
log(1 +7

---

in performance because the number of parameter update for 0 is largely reduced. In the experiment for training epoch number (Figure 4(b)), we can see that the performance is better when the epoch number is larger (the solution of lower-level objective is more precise).
0.90
049
0.91
Best veldattonn
0.8
0.82
K-1
K03
epechass
0.89
Ke5
opoch-
0.75
epoch
0.88
Upper toop iterabonslm
Upper loop iterat onsim!
0.92
0.90
0.92
Det velditon AVE
0.50
0.84
0.89
0.82
008
K-3
spech
Ke5
spoch-
0.87
0.7
200
200
25C
100
250
Running timets
(a) number sample
(b) number of training epochs.
Figure 4: Ablation study for our proposed method on synthetic data.
4.4 Time Cost in Experiments
In this section, evaluate the time cost different parts of our method in three experiments above. The total time cost of these hyper-parameter optimization problems can be mainly divided into two parts: (i) the model training with a specific hyper-parameter and (ii) the update step of 0. Table shows the the clock time in different datasets, where we conduct each experiment for 5 times. As we see, time cost in model training is consistently far more than time cost update step of 0. which shows that it is worthwhile to spend more time to exploit the curvature of upper-level optimization objective in (1).
Table 2: Clock time (in seconds) for model training and update step of 0
synthetic data (Section 4.1)
FB15k237
WNISRR (Section 4.2.1)
CIFAR10 (Section 4.2.2)
(Section 4.2.1)
model training
710 37
86298 + 760
79471 + 936
163992 + 801
update step of 0
16.1 4.2
79.1 +: 2.7
72.4 + 2.3
14.5 + 1.4
Conclusion
In this paper, we address the challenge of hyper- hyper-paramend parameter optimization in scenarios where hyper- gradients are unavailable. To tackle this problem, we introduce novel algorithm that leverages cubic regularization based on stochastic relaxation. Our proposed method offers several advantages over existing approaches, including accelerated convergence and enhanced capability to escape from undesirable saddle points. Extensive experiments conducted both synthetic and real-world datasets demonstrate the effectiveness of our proposed method in solving wide range of hyper- pareeeer optimization problems.

---

Acknowledgments
This work was supported in part by The National Key Research and Development Program of China under grant 2020AA grant from the Research Grants Council of the Hong Kong Special Administrative Region, China (Project No. HKU C7004-22G), and NSF of China (No. 92270106).
References
[1] Ron Kohavi. A study of cross-validation and bootstrap for accuracy estimation and model selection. In International Joint Conference on Artificial Intelligence, 1995.
[2] Jasper Snoek, Hugo Larochelle, and Ryan P Adams. Practical bayesian optimization of machine learning algorithms. Advances in Neural Information Processing Systems, 2012.
[3] James Bergstra and Yoshua Bengio. Random search for of Machine Learning Research, 13(2):281-305 2012.
optimization. Journal
[4] Luca Franceschi, Paolo Frasconi, Saverio Salzo, Riccardo Grazzi, and Massimiliano Pontil. Bilevel programming for hyperparameter optimization and meta-learning. In International Conference Machine Learning, 2018.
[5] Riccardo Grazzi, Luca Franceschi, Massimiliano Pontil, and Saverio Salzo, On the iteration complexity of hypergradient computation. In International Conference on Machine Learning. 2020.
[6] Valerii Likhosherstov, Xingyou Song, Krzysztof Choromanski, Jared Q Davis, and Adrian Weller. Debiasing first-order heuristic for approximate bi-level optimization. In International Conference on Machine Learning, 2021.
[7] Kaiyi Ji, Junjie Yang, and Yingbin Liang. Bilevel optimization: Convergence analysis and enhanced design. In International Conference on Machine Learning, 2021.
[8] Théo Trouillon, Johannes Welbl, Sebastian Riedel, Éric Gaussier, and Guillaume Bouchard. Complex embeddings for simple link prediction. In International Conference on Machine Learning, 2016.
[9] Quan Wang, Zhendong Mao, Bin Wang, and Li Guo. Knowledge graph embedding: A survey of approaches and applications. IEEE Transactions on Knowledge and Data Engineering. 29(12):2724-2744 2017.
[10] Hong-Jian Xue, Xinyu Dai, Jianbing Zhang, Shujian Huang, and Jiajun Chen. Deep matrix factorization models for recommender systems. In International Joint Conference on Artificial Intelligence, 2017.
[11] James Bergstra, Rémi Bardenet, Yoshua Bengio, and Balázs Kégl. Algorithms for hyper- parameter optimization. Advances in Neural Information Processing Systems, 2011.
[12] Lingxi Xie and Alan Yuille. Genetic CNN. In IEEE International Conference on Computer Vision, 2017.
[13] Stuart Geman and Donald Geman. Stochastic relaxation, gibbs distributions, and the bayesian restoration of images. IEEE Transactions Pattern Analysis and Machine Intelligence, 6(6):721-741, 1984.
[14] Sirui Xie, Hehui Zheng, Chunxiao Liu, and Liang Lin. Snas: stochastic neural architecture search. In International Conference on Learning Representations, 2019.
[15] Youhei Akimoto, Shinichi Shirakawa, Nozomu Yoshinari, Kento Uchida, Shota Saito, and Kouhei Nishida. Adaptive stochastic natural gradient method for one-shot neural architecture search. In International Conference on Machine Learning. 2019.
[16) Barret Zoph and Quoe V. Le. Neural architecture search with reinforcement learning. In International Conference on Learning Representations, 2017.
10

---

[17] Yurii Nesterov and Boris T Polyak. Cubic regularization of newton method and its global performance. Mathematical Programmingg 108(1):177-205, 2006.
[18] Matthias Feurer and Frank Hutter. Hyperparameter optimization. Aulomated machine learning: Methods, systems, challenges, pages 3-33, 2019.
[19] Li Yang and Abdallah Shami. On hyperparameter optimization of machine learning algorithms: Theory and practice. Neurocomputting, 415:295-316 2020.
[20] Frank Hutter, Lars Kotthoff, and Joaquin Vanschoren, editors. Automated Machine Learning: Methods, Systems, Challenges. Springer, 2018.
[21] Benoft Colson, Patrice Marcotte, and Gilles Savard. An overview of bilevel optimization. Annals of Operations Research, 153(1):235-225)) 2007.
[22] Hugo Larochelle, Dumitru Erhan, Aaron Courville, James Bergstra, and Yoshua Bengio. An empirical evaluation of deep architectures on problems with many factors of variation. In International Conference on Machine Learning. 2007.
[23] Lisha Li, Kevin Jamieson, Giulia DeSalvo, Afshin Rostamizadeh, and Ameet Talwalkar. Hy- perband: A novel bandit-based approach to hyperparameter optimization. Journal of Machine Learning Research, 18(1):6765-681 161,, 2017.
[24] Coralia Cartis, Nicholas IM Gould, and Philippe L Toint. Adaptive cubic regularisation methods for unconstrained optimization, part i: motivation, convergence and numerical results. Mathematical Programming, 127(2):245-299)) 2011.
[25] Jonas Moritz Kohler and Aurelien Lucchi. Sub-sampled cubic regularization for non-convex optimization. In International Conference on Machine Learning, 2017.
[26] Nilesh Tripuraneni, Mitchell Stern, Chi Jin, Jeffrey Regier, and Michael Jordan. Stochastic cu- bic regularization for fast nonconvex optimization. Advances in Neural Information Processing Systems, 2018
[27] Zhe Wang. Yi Zhou, Yingbin Liang, and Guanghui Lan. Stochastic variance-reduced cubic regularization for nonconvex optimization. In International Conference on Artificial Intelligence and Statistics, 2019.
[28] Yongqi Zhang, Zhanke Zhou, Quanming Yao, and Yong Li. Kgtuner: Efficient hyper-parameter search for knowledge graph learning, arXiv preprint arXiv:2205.02460 2022.
[29] Bo Han, Quanming Yao, Xingrui Yu, Gang Niu, Miao Xu, Weihua Hu, Ivor Tsang, and Masashi Sugiyama. Co-teaching: Robust training of deep neural networks with extremely noisy labels. Advances in Neural Information Processing Systems, 2018.
[30] Quanming Yao, Hansi Yang, Bo Han, Gang Niu, and James Kwok. Searching to exploit memorization effect in learning from corrupted labels. In International Conference on Machine Learning, 2020. 2020.
[31] Minhui Huang, Kaiyi Ji, Shiqian Ma, and Lifeng Lai. Efficiently escaping saddle points in bilevel optimization. arXiv preprint arXiv:2202.03664 2022.
[32] Isabelle Guyon and André Elisseeff. An introduction to variable and feature selection. Journal of machine learning research, 315 2003.
[33] Yutaro Yamada, Ofir Lindenbaum, Sahand Negahban, and Yuval Kluger. Feature selection using stochastic gates. International Conference on Machine Learning, 2020.
[34] Hanxiao Liu, Karen Simonyan, and Yiming Yang. Darts: Differentiable architecture search. In International Conference Learning Representations, 2019.
[35] Shun-Ichi Amari. Natural gradient works efficiently in learning. Neural computation, 10(2):251-
276, 1998.

---

[36] Peng Xu, Fred Roosta, and Michael W Mahoney. Newton-type methods for non-convex optimization under inexact hessian information. Mathematical Programming. 184(1-2):35-70, 2020.
[37] Hadi S. Jomaa, Josif Grabocka, and Lars Schmidt-Thieme. Hyp-rl: Hyperparameter optimiza- tion by reinforcement learning. arXiv preprint 11XXv:/906.115.... 2019.
[38] Kristina Toutanova and Danqi Chen. Observed versus latent features for knowledge base and text inference. In Proceedings of the 3rd Workshop Continuous Vector Space Models and Their Compositionaliit, 2015.
[39] Tim Dettmers, Pasquale Minervini, Pontus Stenetorp, and Sebastian Riedel. Convolutional 2d knowledge graph embeddings. In AAAI Conference on Artificial Intelligence, 2018.
[40] Joel A. Tropp. An introduction to matrix concentration inequalities. Foundations and Trends in Machine Learning, 8(1-2)11-230, 2015.
[41] Prashant Khanduri, Siliang Zeng, Mingyi Hong, Hoi-To Wai, Zhaoran Wang, and Zhuoran Yang. near-optimal algorithm for stochastic bilevel optimization via double- momentummm Advances in Neural Information Processing Systems, 2021.
[42] Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv: 2014.
[43] Bishan Yang, Scott Wen-tau Yih, Xiaodong He, Jianfeng Gao, and Li Deng. Embedding entities and relations for learning and inference in knowledge bases. In International Conference on Learning Representations, 2015.
[44] Timothée Lacroix, Nicolas Usunier, and Guillaume Obozinski. Canonical tensor decomposition for knowledge base completion. In International Conference on Machine Learning, 2018.
[45] Zhanqiu Zhang, Jianyu Cai, and Jie Wang. -uality-iduueed regularizer for tensor factorization based knowledge graph completion. Advances in Neural Information Processing Systems, 2020.
[46] Antoine Bordes, Nicolas Usunier, Alberto Garcia-Duran, Jason Weston, and Oksana Yakhnenko. Translating embeddings for modeling multi-relational data. Advances in Neural Information Processing Systems, 2013.
[47] Zhiqing Sun, Zhi- Zhi-Hong Deng, Jian- Yun Nie, and Jian Tang. Rotate: Knowledge graph em- bedding by relational rotation in complex space. In International Conference on Learning Representations, 2019.
[48] Théo Trouillon, Christopher R Dance, Éric Gaussier, Johannes Welbl, Sebastian Riedel, and Guillaume Bouchard. Knowledge graph completion via complex tensor factorization. Journal of Machine Learning Research, 18: 1-:1-38, 2017.
[49] Zhiqing Sun, Shikhar Vashishth, Soumya Sanyal, Partha Talukdar, and Yiming Yang. A re- evaluation of knowledge graph completion methods. Annual Meeting of the Association for Computational Linguistics, 2020.
12


Reviews:
[{"review_id": "WVBUeONvUC", "reviewer": "Reviewer_S7HE", "summary": "This paper proposes a new optimization based technique for hyper-parameter tunings using Adaptive Cubic Regularized Newton method (ARC) based on stochastic relaxation. The author highlights the limitations of the existing hyper-parameter optimization algorithms.They show that their suggested method achieves better convergence guarantees as existing methods only prove convergence to first order methods where their work also achieves second order optimality conditions. The authors also run some numerical experiments to show that the algorithm outperforms existing methods in terms of faster convergence. ", "questions": "* A typo under Table 1 should be fixed.\n* In all figures, the x-axis is confusing. What do you mean by number of trained models? Also does this mean that the AUC on the y-axis is computed by averaging across all the trained models?\n* What machine learning model you used for all the experiments. More details are needed for example for the classification did you use Decision Tree, SVM, Logistic regression …. Or in the image setup, what neural network models you used?\n* In figure 4, you can add (M) next to the upper loop iterations so that it is easier to understand the plot and what parameters are you considering in your ablation study. Also add the unit for the relative training time.\n* The details about the experiment are important as reproducibility of the work is a major concern. Even if you don’t add a new dataset or create a new machine learning model, there is still a need to reproduce your results and validate that. \n* The details of the experiments should include: what ML models are being trained, the environment, how to set up the datasets and access them, which packages are being used, ARC solver package ….\n* The code of your work also should be attached to the submission and later in your final paper, it is recommended to add the github repository link.\n* Also in your code, did you fix the seed so that people can get similar results as what you got.\n\nComments about the proof:\n* I checked the proof carefully and most of the math seems correct, however following the proof is a bit tricky and I would suggest it if you can make it more clear. Also, if you can give more explanation about the intuition behind those lemmas. \n* Some of the lemmas better to be referenced as facts since you are using the results of previous works: Lemma B.1, Lemma B.2, Lemma B.4. \n* In Lemma B.1 if you can refer to the Hessian as $\\rho$-Lipschitz.\n* Prerequisite of Theorem 3.6 should appear on its own as an equation because you are referring to it in multiple places.\n* In the proof of Preposition 3.2, there is a typo. There should be $\\grad \\log$ in the last two lines of the proof.\n* In Lemma B.3, if you can elaborate more the steps as it is not clear how you get rid of the summation term that appears in the definition of the approximate gradient $g^m$. Since the expectation is with respect to the probability $p_\\theta$ and your summation is in term of $1/k$, but the $p_\\theta$ isn’t uniform so to interchange those terms.\n* In Lemma B.3, how did you get rid of the $\\min$ term that appears in the prerequisite of Theorem 3.6 since you only used the first term in from $\\min (term_1, term_2)$.\n* In Lemma B.5, similar comment to Lemma B.3. How you get rid of the summation term that appears in the definition of the approximate hessian. Also, you only used the second term that appears in the $\\min(term_1, term_2)$ from the prerequisite of Theorem 3.6. \n* In Lemma B.5, in your first line of proof you say “By definition, we have“, but you didn’t refer to any definition. \n* In Lemma B.5, I guess you have a typo in the equality that involves $K$. I guess, this should be $b_2 ^ 2$ instead of $c_2 ^ 2$.\n* In Lemma B.5, I guess you have a typo in the inequality bounding the norm of the difference between the approximate hessian and the actual hessian. $\\|| B ^ m - \\grad J(\\theta ^ m)\\||$ should be  $\\|| B ^ m - \\grad ^ 2 J(\\theta ^ m)\\||$. Also in this lemma, you use the l-2 norm and the spectral norm interchangeably but it is not clear if the math is adjusted at all steps. Or you have some typos.\n* Lemma B.6 is well known facts about ARC, so you can refer to some citations and in that case there is no need to do the proof.\n* In Lemma B.7, when you refer to both Lemmas B.3 and B.5 you mention the inequalities from those lemmas, but without taking into consideration that each one satisfies the inequality with a different probability.\n* In lemma B.9, similar comment to the above when referring to lemmas B.3 and B.5.\n* In lemma B.9, the inequality that appears after the line “And from Lemma B.5, we will have” seems wrong since as I understand you are using Cauchy Schwarz inequality here and it is not clear why you have $\\||\\delta ^ m\\|| ^ 2$ and from where the half term appears. I guess the correct term should be $c_2 \\sqrt{\\rho \\epsilon} \\|| \\delta ^ m \\||$. Then the proof needs to be updated.\n* In Lemma B.9, also it is not clear from where you got the result in case 2 that $\\tilde c(\\delta ^ m)$  <= -\\frac{1}{96} \\sqrt {\\epsilon ^ 3 / \\rho}$.\n", "limitations": "Exposed by the authors.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "* A weakness of existing methods is clearly identified and solved.\n* The originality lies in being the first to use the second order method (ARC) in a bi-level optimization problem for hyper-parameter tuning.\n* The paper is well written and the result can serve as the base of future work in  using second order methods for hyper-parameter tunings.\n* Some minimal numerical results are shown\n", "weaknesses": "* The major contribution of this work should be explicitly stated in a section called “Our contributions”.\n* The authors use ' the curse of dimensionality' as an example of a disadvantage for existing methods; however, when explaining their proposed approach and arguing that even they have to compute the full Hessian, this will not be an expensive operation since hyper-parameters are often low dimensional. This contradicts with the disadvantage mentioned earlier.\n* The authors didn’t compare their results against the global optimal (if exists) set of hyper-parameters in their numerical experiments. \n* In Algorithm.1, the authors didn’t explain how to optimize the lower-level objective in (1)\n* In Algorithm.1, the authors didn’t explain if ARC is being solved using an iterative approach or a factorization based approach when solving the subproblems.\n* In their analysis, the authors are not taking into account the cost of solving the ARC subproblems.\n* For the synthetic data experiments, the authors didn’t use any reference and they didn’t explain why they didn't use real data for the problem of filtering useful features.\n* The details of the experiment parts are not enough so that people can reproduce the results. For example, the authors don't state what machine learning models they used.\n* The proof is a bit tricky and it has some typos and errors.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your comments on both strength and weakness of our paper. Here are our responses to your comments on weakness and questions:\n\n## Weakness\n\n> W1: The major contribution of this work should be explicitly stated in a section called \"Our contributions\".\n\nA1: Please see our response to Q1 in the general response part. \n\n> W2: The authors use 'the curse of dimensionality' as an example of a disadvantage for existing methods; however, when explaining their proposed approach and arguing that even they have to compute the full Hessian, this will not be an expensive operation since hyper-parameters are often low dimensional. This contradicts with the disadvantage mentioned earlier.\n\nA2: Thank you for pointing this out. As the hyper-parameters are often low-dimensional, it may not be appropriate to refer to \"curse of dimension\" for a disadvantage of existing methods. Nevertheless, even with low-dimensional hyper-parameters, existing methods still suffer from slow convergence, as is also demonstrated in experiment part. We will revise this part to avoid such contradictions in final version. \n\n> W3: Compare against global optimal set of hyper-parameters in experiments.\n\nA3: For experiments on real-world data, it is hard to obtain the global optimal set of hyper-parameters in advance, so here we only consider the experiment on synthetic data, where the global optimal set of hyper-parameters will be masking all the noise and retaining all the useful features. We have added its performance to the revised Figure 1(a) and you could check it in the PDF file in general response part. From the revised figure, we can see that the best validation AUC obtained by the proposed method is close to global optimum, while all the other baseline methods have much larger gaps. \n\n> W4: In Algorithm.1, the authors didn’t explain how to optimize the lower-level objective in (1).\n\nA4: Our proposed method should be compatible with any optimizer for lower-level objective, e.g. SGD, Adam or AdamW. In our experiments, we use AdamW in the experiment on synthetic data and use Adam in experiments on real-world data. We will make this more clear in the final version. \n\n> W5: In Algorithm.1, the authors didn’t explain if ARC is being solved using an iterative approach or a factorization based approach when solving the subproblems.\n\nA5: In Algorithm 1, ARC is solved using an iterative approach, similar to the method described in section 6 in [1]. We will make this more clear in the final version.\n\n> W6: In their analysis, the authors are not taking into account the cost of solving the ARC subproblems.\n\nA6: Thank you for raising this question. In our experiments, the number of hyper-parameters is often much less than 100, while the models in lower-level problems have millions of parameters. As such, it should take much longer time to solve lower-level problems than solve the ARC subproblems. It is also empiricaly verified in Section 4.4 that the time cost of updating $\\theta$ (which includes the time cost of solving the ARC subproblems) is negligible in most cases. Therefore, we focus more on the time cost on solving lower-level problems, rather than solving the ARC subproblems. \n\n> W7: For the synthetic data experiments, the authors didn’t use any reference and they didn’t explain why they didn't use real data for the problem of filtering useful features.\n\nA7: We will add citations on feature selection problem in final version, e.g. [2], [3]. As we already have experiment for real data in Section 4.2, the use of synthetic data also enables us to compare against global optimal, which may be hard to obtain if we use real-world data. \n\n> W8: Experiment details\n\nA8: Please see our response to Q3 in common responses.\n\n> W9: Typos and errors in proof\n\nA9: We will thoroughly revise our paper and appendix in our final version. \n\n[1] Adaptive cubic regularisation methods for unconstrained optimization. Part I: motivation, convergence and numerical results. Mathematical Programming, 2011.\n\n[2] An introduction to variable and feature selection. JMLR, 2003.\n\n[3] Feature selection using stochastic gates. ICML 2020\n\n## Questions\n\n> Q10: Typo under Table 1\n\nA10: Thank you for pointing this out. We will revise it in the final version.\n\n> Q11: In all figures, the x-axis is confusing. What do you mean by number of trained models? Also does this mean that the AUC on the y-axis is computed by averaging across all the trained models?\n\nA11: Thank you for pointing this out. We have changed the x-axis to running time and you could check the revised figures in the PDF file in general responses. The y-axis plot the best validation performances obtained at that time. We will make this more clear in next version. \n\n> Q12: Experiment details\n\nA12: Please see our response to Q3 in common responses.\n\n> Q13: In figure 4, you can add (M) next to the upper loop iterations so that it is easier to understand the plot and what parameters are you considering in your ablation study. Also add the unit for the relative training time.\n\nA13: Thank you for your suggestions. We have uploaded revised figures on running time in the PDF file in general response. We will also add notations to figures on upper loop iterations.\n\n> Q14: The code of your work also should be attached to the submission and later in your final paper, it is recommended to add the github repository link.\n\nA14: Thank you for your suggestion. The code still needs some clean up to make it easier to use. We will make our code public on GitHub along with our final version. \n\n> Q15: Random seed in our code\n\nA15: We have fixed the seeds in all our experiments. We will explicitly mention this in the next version of our paper. \n\n> Q16: Comments about the proof. \n\nA16: Thank you for your thorough comments. We have checked our proof and found that these typos/errors did not affect the correctness of key theoretical results. We will make all these revisions in final version. \n"}, {"review_id": "x8s4Tqp1H6", "reviewer": "Reviewer_gbiM", "summary": "This paper studies the problem of hyper-parameter optimization in the context of machine learning. Through stochastic relaxation, the problem can be formulated as bilevel optimization, and the authors propose to use cubic regularization to solve the optimization problem. It is shown that under some regularization conditions on the loss function and hyper-parameter distribution, the algorithm efficiently converges to an approximate second-order stationary point. Furthermore, experiments are conducted to verify the effectiveness of the proposed method.", "questions": "1. In Assumption 3.5, the estimation of gradient and Hessian may be biased. How does this affect the final convergence rate?\n2. How does the proposed method compared with other methods that are not based on stochastic relaxation?", "limitations": "The authors have adequately addressed the limitations and potential societal impact of their work.", "rating": 7, "confidence": 2, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "1. The paper is well-written and quite easy to read. The authors explain the motivations, method and results in a very clear way.\n2. Extensive experiments are conducted to demonstrate the wide applicability of the proposed method. All experiment settings are explained in details.", "weaknesses": "While the authors provide theoretical guarantees for the proposed method, it seems that these guarantees can be directly derived from known results on the convergence of cubic regularization methods. At the same time, because the problem considered is bilevel, some assumption seems less natural. For example, the Lipschitz-Hessian assumption on $J(\\theta)$ is quite difficult to verify in practice, since the expression of $J$ itself contains a minimization problem. I wonder if it is possible to obtain convergence results under more 'fundamental' assumptions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your comments for both the strength and weakness on our paper. We would like to first answer the questions raised in the weakness part: \n\n> W1: While the authors provide theoretical guarantees for the proposed method, it seems that these guarantees can be directly derived from known results on the convergence of cubic regularization methods. \n\nA1: It is not trivial to obtain our theoretical analysis from known results on the convergence of cubic regularization methods. \nSince we are solving a bi-level optimization problem, the gradient and Hessian used in the cubic regularization methods may be inexact. As such, we need to bound the error of estimated gradient and Hessian from the solution to lower-level problem, which is different from any known results to the best of our knowledge. We will make the technical novelty more clear in our next version. \n\n> W2: some assumption seems less natural. For example, the Lipschitz-Hessian assumption on $J(\\theta)$ is quite difficult to verify in practice, since the expression of $J$ itself contains a minimization problem. I wonder if it is possible to obtain convergence results under more 'fundamental' assumptions.\n\nA2: Thank you for your suggestion. Recall that the definition of Lipshitz Hessian is $|| \\nabla^2 J(\\theta) - \\nabla^2 J(\\phi) || \\le \\rho || \\theta -\\phi ||$ for any $\\theta$ and $\\phi$, and from Proposition 3.3, we have \n$\\nabla^2 J(\\theta)  =  \\int H^*(\\theta, z) p_{\\theta}(z) dz$, \nwhere $H^*(\\theta, z)$ depends on the objective value, $\\nabla \\log p_{\\theta}(z)$ and $\\nabla^2 \\log p_{\\theta}(z)$. \nThen the assumption of Lipschitz Hessian can be deduced by assuming 1) bounded objective value; and 2) bounded differences of $\\nabla \\log p_{\\theta}(z)$ and $\\nabla^2 \\log p_{\\theta}(z)$ with different $\\theta$, where similar assumptions are made in Assumption 3.5(ii). We will consider changing to more \"fundamental\" assumptions in the next version. \n\n### Answers to remaining questions\n\n> Q3: In Assumption 3.5, the estimation of gradient and Hessian may be biased. How does this affect the final convergence rate?\n\nA3: In our proof, the estimation errors of gradient and Hessian do affect the convergence rate, and inaccurate gradient or Hessian slow down convergence. Nevertheless, the dependency is very complex, and we consider the worst case performance in theorem 3.6 to make it simple.\n\n> Q4: How does the proposed method compared with other methods that are not based on stochastic relaxation?\n\nA4: In Section 4.2 (Experiments on Real-world Data), we compare the proposed method with other methods that are not based on stochastic relaxation, including random search, Bayesian optimization, Hyperband and reinforcement learning. In the two applications, the proposed method achieves the best performance. \n"}, {"review_id": "H8voYG5t5L", "reviewer": "Reviewer_EwuD", "summary": "The authors propose a cubic regularization scheme for the outer loop optimization of bilevel optimization problems. The cubic regularization is very appealing since it can \"avoid saddle points\". Theoretically, authors extend the convergence results of the cubic regularization to inexact gradients to make it work for bilevel optimization. Empirically, authors show that their method achieves better performance on bilevel optimization tasks. Interestingly, they investigate the eigenvalues of the hessian.\n", "questions": "\n- \"These methods [Bayesian optimization, Reinforcement learning] can easily get trapped in saddle points with poor performance\" Could you provide references for this statement? I am not sure to understand why a Bayesian optimization algorithm would be stuck in a saddle point. From what I understand from the paper, hyper-gradient methods with vanilla first-order methods can be trapped in saddle point, and cubic regularization is here to help.\n- What is the point of Lemma 3.1?\n- Definition of what is called a saddle point could be recalled (i.e. negative eigenvalues of the hessian?)\n- How does the proposed method compare to the use of l-BFGS once the hyper-gradient has been computed? (as done in [1] for instance)\n- Experimental results. In my experience, these kinds of results are usually very sensitive to the number of inner steps to solve the optimization problem. In particular, if the estimation of the inner problem is crude, the value of the hyper-gradient can be poorly estimated. Moreover, to my knowledge,  quasi-second-order methods, combined with approximated gradients, can yield poor results in practice. I did not see this crucial hyper-hyper-parameter reported.\n\n- Figure 1: could you detail how exactly this \"relative eigenvalue\" metric is computed? It seems crucial to make the point of the authors but it is not detailed. I guess this is the ratio of the eigenvalues over the maximal eigenvalue.\n\n- Section 4. Could you write exactly the maths of each bilevel problem which is solved? At least in the appendix?\n\n- I would love to see if the proposed algorithm can improve out-of-distribution generalization in meta-learning settings, like [2])\n\n[1] Deledalle, C.A., Vaiter, S., Fadili, J. and Peyré, G., 2014. Stein Unbiased GrAdient estimator of the Risk (SUGAR) for multiple parameter selection. SIAM Journal on Imaging Sciences,\n\n[2] Lee, K., Maji, S., Ravichandran, A. and Soatto, S., 2019. Meta-learning with differentiable convex optimization. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition (pp. 10657-10665).\n\nTypos:\n- Figure 1 eigen value >> eigenvalue\n", "limitations": "NA", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The idea is very interesting and tackles the very important and hard problem of finding a proper outer procedure for bilevel optimization. Avoiding saddle points seems to be a very appealing property!", "weaknesses": "However, it feels like the writing of the paper could be polished (there are unmerged commit line 264 in the supplementary material.). A lot of statements are not properly backed. I am not sure to properly understand the interest of the proposed method.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your comments that our idea is interesting. Here are our responses to the weakness and questions. \n\n### Weakness\n\n> W1: A lot of statements are not properly backed. \n\nA1: The proposed method is backed by (i) theoretical analysis, which demonstrate that it can converge faster than existing gradient-based methods and converge to second-order approximate stationary points; and (ii) empirical results, where experiments on both synthetic and real data sets demonstrate its superior performance compared to existing hyper-parameter optimization methods. We are also willing to provide more support for any statements in our paper. \n\n> W2: I am not sure to properly understand the interest of the proposed method.\n\nA2: The proposed method mainly targets on hyper-parameter optimization problems where the hyper-parameters are discrete or performance metrics are non-differentiable. Please also see our response to Q1 in common responses for the major contributions of our work.\n\n### Questions\n\n> Q3: References for statement: \"These methods [Bayesian optimization, Reinforcement learning] can easily get trapped in saddle points with poor performance\" \n\nA3: Thank you for pointing this out. The problem of stuck in saddle points should be only for gradient-based reinforcement learning. We will update this part in the final version to avoid such ambiguity. \n\n> Q4: What is the point of Lemma 3.1?\n\nA4: Initially, the upper-level objective of hyper-parameter optimization problems is to find the best hyper-parameter $z$ (i.e. optimize $z$). In Eq (1), we transform the upper level objective to optimize $\\theta$ instead of $z$. The point of Lemma 3.1 is to justify that such transformation is reasonable, as optimizing $\\theta$ is equivalent to optimizing $z$. We will make this point more clear in the final version. \n\n> Q5: Definition of what is called a saddle point could be recalled (i.e. negative eigenvalues of the hessian?)\n\nA5: Thank you for mentioning this. We will add this definition next to the definition for second-order stationary point in the final version. \n\n> Q6:  Compare to l-BFGS\n\nA6: We have added experiments using L-BFGS and include it in the revised figure. Please check the PDF file in general response. The performance of L-BFGS is similar to Newton's method in our experiments. We will also add the revised figures to the final version.\n\n> Q7: Experimental results. In my experience, these kinds of results are usually very sensitive to the number of inner steps to solve the optimization problem. In particular, if the estimation of the inner problem is crude, the value of the hyper-gradient can be poorly estimated. Moreover, to my knowledge, quasi-second-order methods, combined with approximated gradients, can yield poor results in practice. I did not see this crucial hyper-hyper-parameter reported.\n\nA7: For experiment on synthetic data, the number of inner steps is set to 40 by default, and we compare the impact of using different number of inner steps in Section 4.3. Results demonstrate that only using few steps in inner loop yields poorly estimated hyper-gradients and worse final performance. When the number of inner steps are too large, it leads to huge computation cost. We will add more descriptions about the number of inner steps in real-world experiments in final version. \n\n> Q8: Figure 1: could you detail how exactly this \"relative eigenvalue\" metric is computed? It seems crucial to make the point of the authors but it is not detailed. I guess this is the ratio of the eigenvalues over the maximal eigenvalue.\n\nA8: Thank you for pointing this out. Your understanding is correct. \"Relative eigenvalue\" refers to the ratio of eigenvalues over the max eigenvalue. We will add more descriptions about it in the final version. \n\n> Q9: Section 4. Could you write exactly the maths of each bilevel problem which is solved? At least in the appendix?\n\nA9: Here we give the mathematical formulation for each problem in our experiments: \n\n(1) For experiment on synthetic data, the mathematical formulation can be referred to Eq (1) in Section 3.1, where $z$ corresponds to feature mask,  $p_{\\theta}$ uses the sigmoid function and details can be found in Section 4.1. We use AUC as the validation evaluation metric. \n\n(2) For hyper-parameter optimization for knowledge graph completion, the mathematical formulation can also be referred to Eq (1), where $z$ corresponds to a set of categorical hyper-parameters (details in Table 1), $p_\\{\\theta\\}$ uses the softmax-like function as defined in Section 4.2.1. We use MRR as the validation evaluation metric. \n\n(3) For schedule search on learning with noisy training data, the mathematical formulation can be writen as: $\\theta^*=\\arg \\min_{\\theta} J(\\theta), s.t. w^*(z)=\\arg \\min_{w} L_{tra}(w,R_z)$, where $J(\\theta)=E_{z \\sim p_{\\theta}(z)} M_{val}(w^*(R_z))$. $R_z$ is the schedule function parameterized by Eq (7) in Section 4.2.2, and $z$ is the corresponding hyperparameter. For $p_\\{\\theta\\}$, we use Beta distribution for each element of $\\beta$ in Eq (7) and use Dirichlet distribution for $\\alpha$ in Eq (7). The validation evaluation metric is accuracy.\n\nWe will add these formulations for each problem to Appendix in final version. \n\n> Q10: I would love to see if the proposed algorithm can improve out-of-distribution generalization in meta-learning settings, like [2]\n\nA10: Thank you for raising this interesting application. As far as we know, most meta-learning problems deal with differentiable evaluation metrics and continuous hyper-parameters (e.g., model initialization). Therefore, they are not the main focus of our method. We will see if there are some problems in out-of-distribution generalization with non-differentiable metrics or discrete hyper-parameters, where we believe the proposed method should also lead to improvements. \n\n> Q11: Writing and typos. \n\nA11: Thank you for mentioning these errors. We will revise all of them in the final version of the paper. "}, {"review_id": "AP3S0c4fWs", "reviewer": "Reviewer_yqty", "summary": "{Summary}\nThis paper proposes a stochastic cubic regularization type algorithm for hyperparameter optimization that does not depend on hyper-gradients. Theoretical analysis shows that the proposed method can converge to approximate second order stationary points with lower sample complexity than that of first-order optimization methods, which can only find first-order stationary points. Experiments demonstrate the effectiveness of the proposed method using both synthetic and real data.", "questions": "see the previous section.", "limitations": "see the previous section.", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 2, "strengths": "{Strengths}\n1. This paper is clearly written and provides a new algorithm to solve nonconvex bilevel hyperparameter optimization problems (with stochastic relaxation).  \n2. In this setting, solving the cubic subproblem takes little time, as the cubic problem dimension is usually very small.\n3. Theoretical analysis shows that the proposed cubic algorithm achieves lower sample complexity than that of first-order methods.\n\n", "weaknesses": "{Weakness and Questions}\n1. The algorithm design is a direct application of cubic regularization, and therefore is not novel.\n2. It is not clear what is the technical novelty in the convergence analysis compare to the analysis of inexact cubic regularization. The author mentions that the constructed $g$ and $B$ are not unbiased estimators of the gradient and Hessian. How does this challenge and affect the technical proof?\n3. There are existing works on finding second-order stationary points of bi-level optimization, e.g., ``Efficiently Escaping Saddle Points in Bilevel Optimization''. Please discuss and compare with it.\n4.  In the experiments, many curves are piece-wise flat (constant). Does that mean the hyperparameters do not change in those iterations?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your comments on both the strength and weakness of our paper. For your questions raised in the weakness part, here are our answers to each of them:\n\n> Q1: The algorithm design is a direct application of cubic regularization, and therefore is not novel.\n\nA1: Our algorithm is different from existing works from the following perspectives: \n1) in implementation, we need to compute gradient and Hessian from solving the lower-level problem, as in step 8 in Algorithm 1; \n2) theory, we consider optimization with inexact gradient and Hessian, which is also different from directly using cubic regularization.\n \n> Q2: It is not clear what is the technical novelty in the convergence analysis compare to the analysis of inexact cubic regularization. \n\nA2: Thank you for raising this question. Our convergence analysis is different from existing analysis for inexact cubic regularization methods as the inexactness here comes from solving the lower-level problem. As such, we need to bound the error of estimated gradient and Hessian from the solution to lower-level problem, while existing analysis directly assumes error bound on estimating gradient or Hessian. \n\n> Q3: The author mentions that the constructed g and B are not unbiased estimators of the gradient and Hessian. How does this challenge and affect the technical proof?\n\nA3: As g (resp. B) are not unbiased estimators for gradient (resp. Hessian), the estimation error now contains two parts: stochastic noise and inexactness in solving lower-level problems. These parts need to be separately bounded in our proof. We will mention these differences in the final version of our paper. \n\n> Q4: There are existing works on finding second-order stationary points of bi-level optimization, e.g., ``Efficiently Escaping Saddle Points in Bilevel Optimization''. Please discuss and compare with it.\n\nA4: Thank you for pointing this out. We are aware of this paper and it is based on hyper-gradient for solving bi-level optimization problems. As such, it cannot be directly applied to our settings here as the hyper-gradient is unavailable due to discrete hyper-parameters or non-differentiable evaluation metrics (e.g., MRR in knowledge completion task, see Section 4.2.1). We will add some discussions on this paper in our final version. \n\n> Q5: In the experiments, many curves are piece-wise flat (constant). Does that mean the hyperparameters do not change in those iterations?\n\nA5: In our experiments, we report the best validation performances of found hyper-parameters. (similar figures can be referred to Figure 1 in [1] or Figure 2 in [2]) When curves keep piece-wise flat, it means that model with newly sampled hyper-parameter does not obtain better performance than the current optimum. We will add more descriptions to avoid any misunderstandings. \n\n[1] Initializing bayesian hyperparameter optimization via meta-learning, Proceedings of the AAAI Conference on Artificial Intelligence. \n\n[2] SMAC3: A versatile Bayesian optimization package for hyperparameter optimization, The Journal of Machine Learning Research.\n"}]

Author rebuttal:
We thank all reviewers for their constructive feedback. Here we colloct some common questions and reply them in this general response. A PDF file containing all revised figures is also uploaded. 

> Q1: major contributions

A1: The major contributions of this work include:

(1) We propose to utilize cubic regularization to accelerate convergence and avoid saddle points in hyper-parameter optimization problems.

(2) We provide theoretical analysis of the proposed method, showing that the proposed method converge to approximate second-order stationary points with inexactly solved lower level objective considered. 

(3) We verfied the effectiveness of the proposed method in experiment on synthetic and real-world data. 

We will also update these contributions to a separate "Our contributions" section in the final version of the paper. 

> Q2: typos and errors

A2: Thank you for pointing these out. We will make thorough revisions in final version. 

> Q3: experiment details

A3: Thank you for your suggestions. Most details on our experiments can be found in appendix. Here we present some details for easy reference: 

(1) ML model. For experiment on synthetic data, we use a linear model with single layer as feature classifier. For hyper-parameter optimization for knowledge graph completion, we use traditional embedding based model, e.g., TransE, RotatE, ComplEx (see Table 1). For schedule search on learning with noisy training data, we use a 5-layer CNN similar to LeNet as the model for image classification. We will update more details to revised appendix in final version. 

(2) Environment. For all experiments, we use Python 3.7, PyTorch 1.13 and torchvision 0.14. Installing these packages should be enough to replicate all our experiments. 

(3) Datasets. For synthetic data, the generation process is described in Line 232-236 in Section 4.1: "We construct a synthetic dataset on 5-way classification, where the inputs are 50-dimensional vectors. Of the 50 features, 25 have different distributions based on their ground-truth labels, generated from Gaussian distributions with different means for each class and the same variance. The remaining 25 features are filled with i.i.d. Gaussian white noise." 
For knowledge graph completion, the two data sets, FB15k-237 and WN18RR, are well-known in the field of KG. 
For learning with noisy labels, the CIFAR-10 data set is also well-known and can be easily accessed from torchvision package. 

(4) Packages. Despite standard packages for deep learning (e.g. PyTorch), we do not use extra packages (e.g. ARC solver packages) in our experiment.

We will also update more details to revised appendix in final version. 

> Q4: revise figures

A4: We have uploaded the revised figures in the attached PDF file. Our modifications include: 

(1) Add the performance global optimal in Figure 1(a) for easy reference

(2) Add the result of BFGS in Figure 1 

(3) Adjust descriptions of axes in all figures 

Meta-review: The authors propose using cubic regularization to solve stochastic relaxations of hyperparameter optimization problems. This is not a standard application of cubic regularization because they apply it to a bilevel optimization problem where the lower-level is inexactly solved. Experimental results indicate that the approach is promising.
Decision: Accept (poster)