A Survey of Optimization Methods from a Machine Learning Perspective

I Introduction

Recently, machine learning has grown at a remarkable rate, attracting a great number of researchers and practitioners. It has become one of the most popular research directions and plays a significant role in many fields, such as machine translation, speech recognition, image recognition, recommendation system, etc. Optimization is one of the core components of machine learning. The essence of most machine learning algorithms is to build an optimization model and learn the parameters in the objective function from the given data. In the era of immense data, the effectiveness and efficiency of the numerical optimization algorithms dramatically influence the popularization and application of the machine learning models. In order to promote the development of machine learning, a series of effective optimization methods were put forward, which have improved the performance and efficiency of machine learning methods.

From the perspective of the gradient information in optimization, popular optimization methods can be divided into three categories: first-order optimization methods, which are represented by the widely used stochastic gradient methods; high-order optimization methods, in which Newton’s method is a typical example; and heuristic derivative-free optimization methods, in which the coordinate descent method is a representative.

As the representative of first-order optimization methods, the stochastic gradient descent method , as well as its variants, has been widely used in recent years and is evolving at a high speed. However, many users pay little attention to the characteristics or application scope of these methods. They often adopt them as black box optimizers, which may limit the functionality of the optimization methods. In this paper, we comprehensively introduce the fundamental optimization methods. Particularly, we systematically explain their advantages and disadvantages, their application scope, and the characteristics of their parameters. We hope that the targeted introduction will help users to choose the first-order optimization methods more conveniently and make parameter adjustment more reasonable in the learning process.

Compared with first-order optimization methods, high-order methods converge at a faster speed in which the curvature information makes the search direction more effective. High-order optimizations attract widespread attention but face more challenges. The difficulty in high-order methods lies in the operation and storage of the inverse matrix of the Hessian matrix. To solve this problem, many variants based on Newton’s method have been developed, most of which try to approximate the Hessian matrix through some techniques . In subsequent studies, the stochastic quasi-Newton method and its variants are introduced to extend high-order methods to large-scale data .

Derivative-free optimization methods are mainly used in the case that the derivative of the objective function may not exist or be difficult to calculate. There are two main ideas in derivative-free optimization methods. One is adopting a heuristic search based on empirical rules, and the other is fitting the objective function with samples. Derivative-free optimization methods can also work in conjunction with gradient-based methods.

Most machine learning problems, once formulated, can be solved as optimization problems. Optimization in the fields of deep neural network, reinforcement learning, meta learning, variational inference and Markov chain Monte Carlo encounters different difficulties and challenges. The optimization methods developed in the specific machine learning fields are different, which can be inspiring to the development of general optimization methods.

Deep neural networks (DNNs) have shown great success in pattern recognition and machine learning. There are two very popular NNs, i.e., convolutional neural networks (CNNs) and recurrent neural networks (RNNs), which play important roles in various fields of machine learning. CNNs are feedforward neural networks with convolution calculation. CNNs have been successfully used in many fields such as image processing , video processing and natural language processing (NLP) . RNNs are a kind of sequential model and very active in NLP . Besides, RNNs are also popular in the fields of image processing and video processing . In the field of constrained optimization, RNNs can achieve excellent results . In these works, the parameters of weights in RNNs can be learned by analytical methods, and these methods can find the optimal solution according to the trajectory of the state solution. Stochastic gradient-based algorithms are widely used in deep neural networks . However, various problems are emerging when employing stochastic gradient-based algorithms. For example, the learning rate will be oscillating in the later training stage of some adaptive methods , which may lead to the problem of non-converging. Thus, further optimization algorithms based on variance reduction were proposed to improve the convergence rate . Moreover, combining the stochastic gradient descent and the characteristics of its variants is a possible direction to improve the optimization. Especially, switching an adaptive algorithm to the stochastic gradient descent method can improve the accuracy and convergence speed of the algorithm .

Reinforcement learning (RL) is a branch of machine learning, for which an agent interacts with the environment by trial-and-error mechanism and learns an optimal policy by maximizing cumulative rewards . Deep reinforcement learning combines the RL and deep learning techniques, and enables the RL agent to have a good perception of its environment. Recent research has shown that deep learning can be applied to learn a useful representation for reinforcement learning problems . Stochastic optimization algorithms are commonly used in RL and deep RL models.

Meta learning has recently become very popular in the field of machine learning. The goal of meta learning is to design a model that can efficiently adapt to the new environment with as few samples as possible. The application of meta learning in supervised learning can solve the few-shot learning problems . In general, meta learning methods can be summarized into the following three types : metric-based methods , model-based methods and optimization-based methods . We will describe the details of optimization-based meta learning methods in the subsequent sections.

Variational inference is a useful approximation method which aims to approximate the posterior distributions in Bayesian machine learning. It can be considered as an optimization problem. For example, mean-field variational inference uses coordinate ascent to solve this optimization problem . As the amount of data increases continuously, it is not friendly to use the traditional optimization method to handle the variational inference. Thus, the stochastic variational inference was proposed, which introduced natural gradients and extended the variational inference to large-scale data .

Optimization methods have a significative influence on various fields of machine learning. For example, proposed the transformer network using Adam optimization , which is applied to machine translation tasks. proposed super-resolution generative adversarial network for image super resolution, which is also optimized by Adam. proposed Actor-Critic using trust region optimization to solve the deep reinforcement learning on Atari games as well as the MuJoCo environments.

The stochastic optimization method can also be applied to Markov chain Monte Carlo (MCMC) sampling to improve efficiency. In this kind of application, stochastic gradient Hamiltonian Monte Carlo (HMC) is a representative method where the stochastic gradient accelerates the step of gradient update when handling large-scale samples. The noise introduced by the stochastic gradient can be characterized by introducing Gaussian noise and friction terms. Additionally, the deviation caused by HMC discretization can be eliminated by the friction term, and thus the Metropolis-Hasting step can be omitted. The hyper-parameter settings in the HMC will affect the performance of the model. There are some efficient ways to automatically adjust the hyperparameters and improve the performance of the sampler.

The development of optimization brings a lot of contributions to the progress of machine learning. However, there are still many challenges and open problems for optimization problems in machine learning. 1) How to improve optimization performance with insufficient data in deep neural networks is a tricky problem. If there are not enough samples in the training of deep neural networks, it is prone to cause the problem of high variances and overfitting . In addition, non-convex optimization has been one of the difficulties in deep neural networks, which makes the optimization tend to get a locally optimal solution rather than the global optimal solution. 2) For sequential models, the samples are often truncated by batches when the sequence is too long, which will cause deviation. How to analyze the deviation of stochastic optimization in this case and correct it is vital. 3) The stochastic variational inference is graceful and practical, and it is probably a good choice to develop methods of applying high-order gradient information to stochastic variational inference. 4) It may be a great idea to introduce the stochastic technique to the conjugate gradient method to obtain an elegant and powerful optimization algorithm. The detailed techniques to make improvements in the stochastic conjugate gradient is an interesting and challenging problem.

The purpose of this paper is to summarize and analyze classical and modern optimization methods from a machine learning perspective. The remainder of this paper is organized as follows. Section II summarizes the machine learning problems from the perspective of optimization. Section III discusses the classical optimization algorithms and their latest developments in machine learning. Particularly, the recent popular optimization methods including the first and second order optimization algorithms are emphatically introduced. Section IV describes the developments and applications of optimization methods in some specific machine learning fields. Section V presents the challenges and open problems in the optimization methods. Finally, we conclude the whole paper.

II Machine Learning Formulated as Optimization

Almost all machine learning algorithms can be formulated as an optimization problem to find the extremum of an objective function. Building models and constructing reasonable objective functions are the first step in machine learning methods. With the determined objective function, appropriate numerical or analytical optimization methods are usually used to solve the optimization problem.

According to the modeling purpose and the problem to be solved, machine learning algorithms can be divided into supervised learning, semi-supervised learning, unsupervised learning, and reinforcement learning. Particularly, supervised learning is further divided into the classification problem (e.g., sentence classification , image classification , etc.) and regression problem; unsupervised learning is divided into clustering and dimension reduction , among others.

For supervised learning, the goal is to find an optimal mapping function $f(x)$ to minimize the loss function of the training samples,

where $N$ is the number of training samples, $\theta$ is the parameter of the mapping function, $x^{i}$ is the feature vector of the $i$th samples, $y^{i}$ is the corresponding label, and $L$ is the loss function.

There are many kinds of loss functions in supervised learning, such as the square of Euclidean distance, cross-entropy, contrast loss, hinge loss, information gain and so on. For regression problems, the simplest way is using the square of Euclidean distance as the loss function, that is, minimizing square errors on training samples. But the generalization performance of this kind of empirical loss is not necessarily good. Another typical form is structured risk minimization, whose representative method is the support vector machine. On the objective function, regularization items are usually added to alleviate overfitting, e.g., in terms of $L_{2}$ norm,

where $\lambda$ is the compromise parameter, which can be determined through cross-validation.

II-B Optimization Problems in Semi-supervised Learning

Semi-supervised learning (SSL) is the method between supervised and unsupervised learning, which incorporates labeled data and unlabeled data during the training process. It can deal with different tasks including classification tasks , regression tasks , clustering tasks and dimensionality reduction tasks . There are different kinds of semi-supervised learning methods including self-training, generative models, semi-supervised support vector machines (S3VM) , graph-based methods, multi-learning method and others. We take S3VM as an example to introduce the optimization in semi-supervised learning.

S3VM is a learning model that can deal with binary classification problems and only part of the training set in this problem is labeled. Let $D^{l}$ be labeled data which can be represented as $D^{l}=\{\{x^{1},y^{1}\},\{x^{2},y^{2}\},...,\{x^{l},y^{l}\}\}$, and $D^{u}$ be unlabeled data which can be represented as $D^{u}=\{x^{l+1},x^{l+2},...,x^{N}\}$ with $N=l+u$. In order to use the information of unlabeled data, additional constraint on the unlabeled data is added to the original objective of SVM with slack variables $\zeta^{i}$. Specifically, define $\epsilon^{j}$ as the misclassification error of the unlabeled instance if its true label is positive and $z^{j}$ as the misclassification error of the unlabeled instance if its true label is negative. The constraint means to make $\sum_{j=l+1}^{N}\min(\epsilon^{i},\zeta^{i})$ as small as possible. Thus, an S3VM problem can be described as

where $C$ is a penalty coefficient. The optimization problem in S3VM is a mixed-integer problem which is difficult to deal with . There are various methods summarized in to deal with this problem, such as the branch and bound techniques and convex relaxation methods .

II-C Optimization Problems in Unsupervised Learning

Clustering algorithms divide a group of samples into multiple clusters ensuring that the differences between the samples in the same cluster are as small as possible, and samples in different clusters are as different as possible. The optimization problem for the $k$-means clustering algorithm is formulated as minimizing the following loss function:

where $K$ is the number of clusters, ${x}$ is the feature vector of samples, $\mu_{k}$ is the center of cluster $k$, and $S_{k}$ is the sample set of cluster $k$. The implication of this objective function is to make the sum of variances of all clusters as small as possible.

The dimensionality reduction algorithm ensures that the original information from data is retained as much as possible after projecting them into the low-dimensional space. Principal component analysis (PCA) is a typical algorithm of dimensionality reduction methods. The objective of PCA is formulated to minimize the reconstruction error as

where $N$ represents the number of samples, ${x}_{i}$ is a $D$-dimensional vector, ${\overline{x}}^{i}$ is the reconstruction of ${x}^{i}$. ${z}^{i}=\{z_{1}^{i},...,z_{D^{\prime}}^{i}\}$ is the projection of ${x}^{i}$ in $D^{\prime}$-dimensional coordinates. ${e}_{j}$ is the standard orthogonal basis under $D^{\prime}$-dimensional coordinates.

Another common optimization goal in probabilistic models is to find an optimal probability density function of $p({x})$, which maximizes the logarithmic likelihood function (MLE) of the training samples,

In the framework of Bayesian methods, some prior distributions are often assumed on parameter $\theta$, which also has the effect of alleviating overfitting.

II-D Optimization Problems in Reinforcement Learning

Reinforcement learning , unlike supervised learning and unsupervised learning, aims to find an optimal strategy function, whose output varies with the environment. For a deterministic strategy, the mapping function from state $s$ to action $a$ is the learning target. For an uncertain strategy, the probability of executing each action is the learning target. In each state, the action is determined by $a=\pi(s)$, where $\pi(s)$ is the policy function.

The optimization problem in reinforcement learning can be formulated as maximizing the cumulative return after executing a series of actions which are determined by the policy function,

where $V_{\pi}(s)$ is the value function of state $s$ under policy $\pi$, $r$ is the reward, and $\gamma\in$ is the discount factor.

II-E Optimization for Machine Learning

Overall, the main steps of machine learning are to build a model hypothesis, define the objective function, and solve the maximum or minimum of the objective function to determine the parameters of the model. In these three vital steps, the first two steps are the modeling problems of machine learning, and the third step is to solve the desired model by optimization methods.

III Fundamental Optimization Methods and Progresses

From the perspective of gradient information, fundamental optimization methods can be divided into first-order optimization methods, high-order optimization methods and derivative-free optimization methods. These methods have a long history and are constantly evolving. They are progressing in many practical applications and have achieved good performance. Besides these fundamental methods, preconditioning is a useful technique for optimization methods. Applying reasonable preconditioning can reduce the number of iterations and obtain better spectral characteristics. These technologies have been widely used in practice. For the convenience of researchers, we summarize the existing common optimization toolkits in a table at the end of this section.

In the field of machine learning, the most commonly used first-order optimization methods are mainly based on gradient descent. In this section, we introduce some of the representative algorithms along with the development of the gradient descent methods. At the same time, the classical alternating direction method of multipliers and the Frank-Wolfe method in numerical optimization are also introduced.

The gradient descent method is the earliest and most common optimization method. The idea of the gradient descent method is that variables update iteratively in the (opposite) direction of the gradients of the objective function. The update is performed to gradually converge to the optimal value of the objective function. The learning rate $\eta$ determines the step size in each iteration, and thus influences the number of iterations to reach the optimal value .

The steepest descent algorithm is a widely known algorithm. The idea is to select an appropriate search direction in each iteration so that the value of the objective function minimizes the fastest. Gradient descent and steepest descent are not the same, because the direction of the negative gradient does not always descend fastest. Gradient descent is an example of using the Euclidean norm in steepest descent .

Next, we give the formal expression of gradient descent method. For a linear regression model, we assume that $f_{\theta}(x)$ is the function to be learned, $L(\theta)$ is the loss function, and $\theta$ is the parameter to be optimized. The goal is to minimize the loss function with

where $N$ is the number of training samples, $D$ is the number of input features, $x^{i}$ is an independent variable with $x^{i}=(x_{1}^{i},...,x_{D}^{i})$ for $i=1,...,N$ and $y^{i}$ is the target output. The gradient descent alternates the following two steps until it converges:

Derive $L(\theta)$ for $\theta_{j}$ to get the gradient corresponding to each $\theta_{j}$:

Update each $\theta_{j}$ in the negative gradient direction to minimize the risk function:

The gradient descent method is simple to implement. The solution is global optimal when the objective function is convex. It often converges at a slower speed if the variable is closer to the optimal solution, and more careful iterations need to be performed.

In the above linear regression example, note that all the training data are used in each iteration step, so the gradient descent method is also called the batch gradient descent. If the number of samples is $N$ and the dimension of $x$ is $D$, the computation complexity for each iteration will be $O(ND)$. In order to mitigate the cost of computation, some parallelization methods were proposed . However, the cost is still hard to accept when dealing with large-scale data. Thus, the stochastic gradient descent method emerges.

III-A2 Stochastic Gradient Descent

Since the batch gradient descent has high computational complexity in each iteration for large-scale data and does not allow online update, stochastic gradient descent (SGD) was proposed . The idea of stochastic gradient descent is using one sample randomly to update the gradient per iteration, instead of directly calculating the exact value of the gradient. The stochastic gradient is an unbiased estimate of the real gradient . The cost of the stochastic gradient descent algorithm is independent of sample numbers and can achieve sublinear convergence speed . SGD reduces the update time for dealing with large numbers of samples and removes a certain amount of computational redundancy, which significantly accelerates the calculation. In the strong convex problem, SGD can achieve the optimal convergence speed . Meanwhile, it overcomes the disadvantage of batch gradient descent that cannot be used for online learning.

The loss function (8) can be written as the following equation:

If a random sample $i$ is selected in SGD, the loss function will be $L^{*}(\theta)$:

The gradient update in SGD uses the random sample $i$ rather than all samples in each iteration,

Since SGD uses only one sample per iteration, the computation complexity for each iteration is $O(D)$ where $D$ is the number of features. The update rate for each iteration of SGD is much faster than that of batch gradient descent when the number of samples $N$ is large. SGD increases the overall optimization efficiency at the expense of more iterations, but the increased iteration number is insignificant compared with the high computation complexity caused by large numbers of samples. It is possible to use only thousands of samples overall to get the optimal solution even when the sample size is hundreds of thousands. Therefore, compared with batch methods, SGD can effectively reduce the computational complexity and accelerate convergence.

However, one problem in SGD is that the gradient direction oscillates because of additional noise introduced by random selection, and the search process is blind in the solution space. Unlike batch gradient descent which always moves towards the optimal value along the negative direction of the gradient, the variance of gradients in SGD is large and the movement direction in SGD is biased. So, a compromise between the two methods, the mini-batch gradient descent method (MSGD), was proposed .

The MSGD uses $b$ independent identically distributed samples ($b$ is generally in 50 to 256 ) as the sample sets to update the parameters in each iteration. It reduces the variance of the gradients and makes the convergence more stable, which helps to improve the optimization speed. For brevity, we will call MSGD as SGD in the following sections.

As a common feature of stochastic optimization, SGD has a better chance of finding the global optimal solution for complex problems. The deterministic gradient in batch gradient descent may cause the objective function to fall into a local minimum for the multimodal problem. The fluctuation in the SGD helps the objective function jump to another possible minimum. However, the fluctuation in SGD always exists, which may more or less slow down the process of converging.

There are still many details to be noted about the use of SGD in the concrete optimization process , such as the choice of a proper learning rate. A too small learning rate will result in a slower convergence rate, while a too large learning rate will hinder convergence, making loss function fluctuate at the minimum. One way to solve this problem is to set up a predefined list of learning rates or a certain threshold and adjust the learning rate during the learning process . However, these lists or thresholds need to be defined in advance according to the characteristics of the dataset. It is also inappropriate to use the same learning rate for all parameters. If data are sparse and features occur at different frequencies, it is not expected to update the corresponding variables with the same learning rate. A higher learning rate is often expected for less frequently occurring features .

Besides the learning rate, how to avoid the objective function being trapped in infinite numbers of the local minimum is a common challenge. Some work has proved that this difficulty does not come from the local minimum values, but comes from the “saddle point” . The slope of a saddle point is positive in one direction and negative in another direction, and gradient values in all directions are zero. It is an important problem for SGD to escape from these points. Some research about escaping from saddle points were developed .

III-A3 Nesterov Accelerated Gradient Descent

Although SGD is popular and widely used, its learning process is sometimes prolonged. How to adjust the learning rate, how to speed up the convergence, and how to prevent being trapped at a local minimum during the search are worthwhile research directions.

Much work is presented to improve SGD. For example, the momentum idea was proposed to be applied in SGD . The concept of momentum is derived from the mechanics of physics, which simulates the inertia of objects. The idea of applying momentum in SGD is to preserve the influence of the previous update direction on the next iteration to a certain degree. The momentum method can speed up the convergence when dealing with high curvature, small but consistent gradients, or noisy gradients . The momentum algorithm introduces the variable $v$ as the speed, which represents the direction and the rate of the parameter’s movement in the parameter space. The speed is set as the average exponential decay of the negative gradient.

In the gradient descent method, the speed update is $v=\eta\cdot(-\frac{\partial L(\theta)}{\partial(\theta)})$ each time. Using the momentum algorithm, the amount of the update $v$ is not just the amount of gradient descent calculated by $\eta\cdot(-\frac{\partial L(\theta)}{\partial(\theta)})$. It also takes into account the friction factor, which is represented as the previous update $v^{old}$ multiplied by a momentum factor ranging between . Generally, the mass of the object is set to 1. The formulation is expressed as

where $mtm$ is the momentum factor. If the current gradient is parallel to the previous speed $v^{old}$, the previous speed can speed up this search. The proper momentum plays a role in accelerating the convergence when the learning rate is small. If the derivative decays to 0, it will continue to update $v$ to reach equilibrium and will be attenuated by friction. It is beneficial for escaping from the local minimum in the training process so that the search process can converge more quickly . If the current gradient is opposite to the previous update $v^{old}$, the value $v^{old}$ will have a deceleration effect on this search.

The momentum method with a proper momentum factor plays a positive role in reducing the oscillation of convergence when the learning rate is large. How to select the proper size of the momentum factor is also a problem. If the momentum factor is small, it is hard to obtain the effect of improving convergence speed. If the momentum factor is large, the current point may jump out of the optimal value point. Many experiments have empirically verified the most appropriate setting for the momentum factor is 0.9 .

Nesterov Accelerated Gradient Descent (NAG) makes further improvement over the traditional momentum method . In Nesterov momentum, the momentum $v^{old}\cdot mtm$ is added to $\theta$, denoted as $\widetilde{\theta}$. The gradient of $\widetilde{\theta}$ is used when updating. The detailed update formulae for parameters $\theta$ are as follows:

The improvement of Nesterov momentum over momentum is reflected in updating the gradient of the future position instead of the current position. From the update formula, we can find that Nestorov momentum includes more gradient information compared with the traditional momentum method. Note that Nesterov momentum improves the convergence from $O(\frac{1}{k})$ (after $k$ steps) to $O(\frac{1}{k^{2}})$, when not using stochastic optimization .

Another issue worth considering is how to determine the size of the learning rate. It is more likely to occur the oscillation if the search is closer to the optimal point. Thus, the learning rate should be adjusted. The learning rate decay factor $d$ is commonly used in the SGD’s momentum method, which makes the learning rate decrease with the iteration period . The formula of the learning rate decay is defined as

where $\eta_{t}$ is the learning rate at the $t$th iteration, $\eta_{0}$ is the original learning rate, and $d$ is a decimal in $$. As can be seen from the formula, the smaller the$d$is, the slower the decay of the learning rate will be. The learning rate remains unchanged when$d=0$and the learning rate decays fastest when$d=1$.

III-A4 Adaptive Learning Rate Method

The manually regulated learning rate greatly influences the effect of the SGD method. It is a tricky problem for setting an appropriate value of the learning rate . Some adaptive methods were proposed to adjust the learning rate automatically. These methods are free of parameter adjustment, fast to converge, and often achieving not bad results. They are widely used in deep neural networks to deal with optimization problems.

The most straightforward improvement to SGD is AdaGrad . AdaGrad adjusts the learning rate dynamically based on the historical gradient in some previous iterations. The update formulae are as follows:

where $g_{t}$ is the gradient of parameter $\theta$ at iteration $t$, $V_{t}$ is the accumulate historical gradient of parameter $\theta$ at iteration $t$, and ${\theta}_{t}$ is the value of parameter $\theta$ at iteration $t$.

The difference between AdaGrad and gradient descent is that during the parameter update process, the learning rate is no longer fixed, but is computed using all the historical gradients accumulated up to this iteration. One main benefit of AdaGrad is that it eliminates the need to tune the learning rate manually. Most implementations use a default value of 0.01 for $\eta$ in (18).

Although AdaGrad adaptively adjusts the learning rate, it still has two issues. 1) The algorithm still needs to set the global learning rate $\eta$ manually. 2) As the training time increases, the accumulated gradient will become larger and larger, making the learning rate tend to zero, resulting in ineffective parameter update.

AdaGrad was further improved to AdaDelta and RMSProp for solving the problem that the learning rate will eventually go to zero. The idea is to consider not accumulating all historical gradients, but focusing only on the gradients in a window over a period, and using the exponential moving average to calculate the second-order cumulative momentum,

where $\beta$ is the exponential decay parameter. Both RMSProp and AdaDelta have been developed independently around the same time, stemming from the need to resolve the radically diminishing learning rates of AdaGrad.

Adaptive moment estimation (Adam) is another advanced SGD method, which introduces an adaptive learning rate for each parameter. It combines the adaptive learning rate and momentum methods. In addition to storing an exponentially decaying average of past squared gradients $V_{t}$, like AdaDelta and RMSProp, Adam also keeps an exponentially decaying average of past gradients $m_{t}$, similar to the momentum method:

where $\beta_{1}$ and $\beta_{2}$ are exponential decay rates. The final update formula for the parameter $\theta$ is

The default values of $\beta_{1}$, $\beta_{2}$, and $\epsilon$ are suggested to set to 0.9, 0.999, and $10^{-8}$, respectively. Adam works well in practice and compares favorably to other adaptive learning rate algorithms.

III-A5 Variance Reduction Methods

Due to a large amount of redundant information in the training samples, the SGD methods are very popular since they were proposed. However, the stochastic gradient method can only converge at a sublinear rate and the variance of gradient is often very large. How to reduce the variance and improve SGD to the linear convergence has always been an important problem.

Stochastic Average Gradient The stochastic average gradient (SAG) method is a variance reduction method proposed to improve the convergence speed. The SAG algorithm maintains parameter $d$ recording the sum of the $N$ latest gradients $\{g_{i}\}$ in memory where $g_{i}$ is calculated using one sample $i,i\in\{1,...,N\}$. The detailed implementation is to select a sample $i_{t}$ to update $d$ randomly, and use $d$ to update the parameter $\theta$ in iteration $t$:

where the updated item $d$ is calculated by replacing the old gradient $\hat{g}_{i_{t}}$ in $d$ with the new gradient $g_{i_{t}}(\theta_{t-1})$ in iteration $t$, $\alpha$ is a constant representing the learning rate. Thus, each update only needs to calculate the gradient of one sample, not the gradients of all samples. The computational overhead is no different from SGD, but the memory overhead is much larger. This is a typical way of using space for saving time. The SAG has been shown to be a linear convergence algorithm , which is much faster than SGD, and has great advantages over other stochastic gradient algorithms.

However, the SAG method is only applicable to the case where the loss function is smooth and the objective function is convex , such as convex linear prediction problems. In this case, the SAG achieves a faster convergence rate than the SGD. In addition, under some specific problems, it can even deliver better convergence than the standard batch gradient descent.

Stochastic Variance Reduction Gradient Since the SAG method is only applicable to smooth and convex functions and needs to store the gradient of each sample, it is inconvenient to be applied in non-convex neural networks. The stochastic variance reduction gradient (SVRG) method was proposed to improve the performance of optimization in the complex models.

The strategies of SAG and SVRG are related to variance reduction. Compared with SAG, SVRG does not need to maintain all gradients in memory, which means that memory resources are saved, and it can be applied to complex problems efficiently. Experiments have shown that the performance of SVRG is remarkable on a non-convex neural network . There are also many variants of such linear convergence stochastic optimization algorithms, such as the SAGA algorithm .

III-A6 Alternating Direction Method of Multipliers

Augmented Lagrangian multiplier method is a common method to solve optimization problems with linear constraints. Compared with the naive Lagrangian multiplier method, it makes problems easier to solve by adding a penalty term to the objective. Consider the following example,

The augmented Lagrange function for problem (27) is

When solved by the augmented Lagrangian multiplier method, its $t$th step iteration starts from the given $\lambda_{t}$, and the optimization turns out to

Separating the $(x,y)$ sub-problem in (29), the augmented Lagrange multiplier method can be relaxed to the following alternating direction method of multipliers (ADMM) . Its $t$th step iteration starts with the given $(y_{t},\lambda_{t})$, and the details of iterative optimization are as follows:

The penalty parameter $\beta$ has a certain impact on the convergence rate of the ADMM. The larger $\beta$ is, the greater the penalties for the constraint term. In general, a monotonically increasing sequence of $\left\{\beta_{t}\right\}$ can be adopted instead of the fixed $\beta$ . Specifically, an auto-adjustment criterion that automatically adjusts $\left\{\beta_{t}\right\}$ based on the current value of $\left\{x_{t}\right\}$ during the iteration was proposed, and applied for solving some convex optimization problems .

The ADMM method uses the separable operators in the convex optimization problem to divide a large problem into multiple small problems that can be solved in a distributed manner. In theory, the framework of ADMM can solve most of the large-scale optimization problems. However, there are still some problems in practical applications. For example, if we use a stop criterion to determine whether convergence occurs, the original residuals and dual residuals are both related to $\beta$, and $\beta$ with a large value will lead to difficulty in meeting the convergence conditions .

III-A7 Frank-Wolfe Method

In 1956, Frank and Wolfe proposed an algorithm for solving linear constraint problems . The basic idea is to approximate the objective function with a linear function, then solve the linear programming to find the feasible descending direction, and finally make a one-dimensional search along the direction in the feasible domain. This method is also called the approximate linearization method.

Here, we give a simple example of Frank-Wolfe method. Consider the optimization problem,

where $A$ is an $m\times n$ full row rank matrix, and the feasible region is $S=\left\{x|Ax=b,x\geq 0\right\}$. Expand $f(x)$ linearly at $x_{0}$, $f(x)\approx f(x_{0})+\nabla f(x_{0})^{\top}(x-x_{0})$, and substitute it into equation (31). Then we have

Suppose there exist an optimal solution $y_{t}$, and then there must be

So $y_{t}-x_{t}$ is the decreasing direction of $f(x)$ at $x_{t}$. A fetch step of $\lambda_{t}$ updates the search point in a feasible direction. The detailed operation is shown in Algorithm 1.

The algorithm satisfies the following convergence theorem :

(1) $x_{t}$ is the Kuhn-Tucker point of (31) when $\nabla f(x_{t})^{\top}(y_{t}-x_{t})=0$.

(2) Since $y_{t}$ is an optimal solution for problem (33), the vector $d_{t}$ satisfies $d_{t}=y_{t}-x_{t}$ and is the feasible descending direction of $f$ at point $x_{t}$ when $\nabla f(x_{t})^{\top}(y_{t}-x_{t})\neq 0$.

The Frank-Wolfe algorithm is a first-order iterative method for solving convex optimization problems with constrained conditions. It consists of determining the feasible descent direction and calculating the search step size. The algorithm is characterized by fast convergence in early iterations and slower in later phases. When the iterative point is close to the optimal solution, the search direction and the gradient direction of the objective function tend to be orthogonal. Such a direction is not the best downward direction so that the Frank-Wolfe algorithm can be improved and extended in terms of the selection of the descending directions .

III-A8 Summary

We summarize the mentioned first-order optimization methods in terms of properties, advantages, and disadvantages in Table LABEL:3.1.8.

III-B High-Order Methods

The second-order methods can be used for addressing the problem where an objective function is highly non-linear and ill-conditioned. They work effectively by introducing curvature information.

This section begins with introducing the conjugate gradient method, which is a method that only needs first-order derivative information for well-defined quadratic programming, but overcomes the shortcoming of the steepest descent method, and avoids the disadvantages of Newton’s method of storing and calculating the inverse Hessian matrix. But note that when applying it to general optimization problems, the second-order gradient is needed to get an approximation to quadratic programming. Then, the classical quasi-Newton method using second-order information is described. Although the convergence of the algorithm can be guaranteed, the computational process is costly and thus rarely used for solving large machine learning problems. In recent years, with the continuous improvement of high-order optimization methods, more and more high-order methods have been proposed to handle large-scale data by using stochastic techniques . From this perspective, we discuss several high-order methods including the stochastic quasi-Newton method (integrating the second-order information and the stochastic method) and their variants. These algorithms allow us to use high-order methods to process large-scale data.

The conjugate gradient (CG) approach is a very interesting optimization method, which is one of the most effective methods for solving large-scale linear systems of equations. It can also be used for solving nonlinear optimization problems . As we know, the first-order methods are simple but have a slow convergence speed, and the second-order methods need a lot of resources. Conjugate gradient optimization is an intermediate algorithm, which can only utilize the first-order information for some problems but ensures the convergence speed like high-order methods.

Early in the 1960s, a conjugate gradient method for solving a linear system was proposed, which is an alternative to Gaussian elimination . Then in 1964, the conjugate gradient method was extended to handle nonlinear optimization for general functions . For years, many different algorithms have been presented based on this method, some of which have been widely used in practice. The main features of these algorithms are that they have faster convergence speed than steepest descent. Next, we describe the conjugate gradient method.

where $A$ is an $n\times n$ symmetric, positive-definite matrix. The matrix $A$ and vector $b$ are known, and we need to solve the value of $\theta$. The problem (35) can also be considered as an optimization problem that minimizes the quadratic positive definite function,

The above two equations have an identical unique solution. It enables us to regard the conjugate gradient as a method for solving optimization problems.

The gradient of $F(\theta)$ can be obtained by simple calculation, and it equals the residual of the linear system : $r(\theta)=\nabla F(\theta)=A\theta-b.$

Conjugate: Given an $n\times n$ symmetric positive-definite matrix $A$, two non-zero vector $d_{i},d_{j}$ are conjugate with respect to $A$ if

A set of non-zero vector $\{d_{1},d_{2},d_{3},....,d_{n}\}$ is said to be conjugate with respect to $A$ if any two unequal vectors are conjugate with respect to $A$ .

Next, we introduce the detailed derivation of the conjugate gradient method. $\theta_{0}$ is a starting point, $\{d_{t}\}_{t=1}^{n-1}$ is a set of conjugate directions. In general, one can generate the update sequence $\{\theta_{1},\theta_{2},....,\theta_{n}\}$ by a iteration formula:

The step size $\eta_{t}$ can be obtained by a linear search, which means choosing $\eta_{t}$ to minimize the object function $f(\cdot)$ along $\theta_{t}+\eta_{t}d_{t}$. After some calculations (more details in ), the update formula of $\eta_{t}$ is

The search direction $d_{t}$ is obtained by a linear combination of the negative residual and the previous search direction,

where $r_{t}$ can be updated by $r_{t}=r_{t-1}+\eta_{t-1}Ad_{t-1}$. The scalar $\beta_{t}$ is the update parameter, which can be determined by satisfying the requirement that $d_{t}$ and $d_{t-1}$ are conjugate with respect to $A$, i.e., $d_{t}^{\top}Ad_{t-1}=0$. Multiplying both sides of the equation (40) by $d_{t-1}^{\top}A$, one can obtain $\beta_{t}$ by

After several derivations of the above formula according to , the simplified version of $\beta_{t}$ is

The CG method, has a graceful property that generating a new vector $d_{t}$ only using the previous vector $d_{t-1}$, which does not need to know all the previous vectors $d_{0},d_{1},d_{2}\dots d_{t-2}$. The linear conjugate gradient algorithm is shown in Algorithm 2.

III-B2 Quasi-Newton Methods

Gradient descent employs first-order information, but its convergence rate is slow. Thus, the natural idea is to use second-order information, e.g., Newton’s method . The basic idea of Newton’s method is to use both the first-order derivative (gradient) and second-order derivative (Hessian matrix) to approximate the objective function with a quadratic function, and then solve the minimum optimization of the quadratic function. This process is repeated until the updated variable converges.

The one-dimensional Newton’s iteration formula is shown as

where $f$ is the object function. More general, the high-dimensional Newton’s iteration formula is

where $\nabla^{2}f$ is a Hessian matrix of $f$. More precisely, if the learning rate (step size factor) is introduced, the iteration formula is shown as

where $d_{t}$ is the Newton’s direction, $\eta_{t}$ is the step size. This method can be called damping Newton’s method . Geometrically speaking, Newton’s method is to fit the local surface of the current position with a quadratic surface, while the gradient descent method is to fit the current local surface with a plane .

Quasi-Newton Method Newton’s method is an iterative algorithm that requires the computation of the inverse Hessian matrix of the objective function at each step, which makes the storage and computation very expensive. To overcome the expensive storage and computation, an approximate algorithm was considered which is called the quasi-Newton method. The essential idea of the quasi-Newton method is to use a positive definite matrix to approximate the inverse of the Hessian matrix, thus simplifying the complexity of the operation. The quasi-Newton method is one of the most effective methods for solving non-linear optimization problems. Moreover, the second-order gradient is not directly needed in the quasi-Newton method, so it is sometimes more efficient than Newton’s method. In the following section, we will introduce several quasi-Newton methods, in which the Hessian matrix and its inverse matrix are approximated by different methods.

Quasi-Newton Condition We first introduce the quasi-Newton condition. Assuming that the objective function $f$ can be approximated by a quadratic function, we can extend $f(\theta)$ to Taylor series at $\theta=\theta_{t+1}$, i.e.,

Then we can compute the gradient on both sides of the above equation, and obtain

Use $B$ to represent the approximate matrix of the Hessian matrix. Set $s_{t}=\theta_{t+1}-\theta_{t}$, and $u_{t}=\nabla f(\theta_{t+1})-\nabla f(\theta_{t})$. The matrix $B_{t+1}$ is satisfied that

This equation is called the quasi-Newton condition, or secant equation.

The search direction of quasi-Newton method is

where $g_{t}$ is the gradient of $f$, and the update of quasi-Newton is

The step size $\eta_{t}$ is chosen to satisfy the Wolfe conditions, which is a set of inequalities for inexact line searches $\min_{\eta_{t}}f(\theta_{t}+\eta_{t}d_{t})$ . Unlike Newton’s method, quasi-Newton method uses $B_{t}$ to approximate the true Hessian matrix. In the following paragraphs, we will introduce some particular quasi-Newton methods, in which $H_{t}$ is used to express the inverse of $B_{t}$, i.e., $H_{t}=B_{t}^{-1}$.

DFP In the 1950s, a physical scientist, William C. Davidon , proposed a new approach to solve nonlinear problems. Then Fletcher and Powel explained and improved this method, which sparked a lot of research in the late 1960s and early 1970s . DFP is the first quasi-Newton method named after the initials of their three names. The DFP correction formula is one of the most creative inventions in the field of non-linear optimization, shown as below:

BFGS Broyden, Fletcher, Goldfarb and Shanno proposed the BFGS method , in which $B_{t+1}$ is updated according to

The quasi-Newton algorithm still cannot solve large-scale data optimization problem because the method generates a sequence of matrices to approximate the Hessian matrix. Storing these matrices needs to consume computer resources, especially for high-dimensional problems. It is also impossible to retain these matrices in the high-speed storage of computers, restricting its use to even small and midsize problems .

L-BFGS Limited memory quasi-Newton methods, named L-BFGS is an improvement based on the quasi-Newton method, which is feasible in dealing with the high-dimensional situation. The method stores just a few $n$-dimensional vectors, instead of retaining and computing fully dense $n\times n$ approximations of the Hessian . The basic idea of L-BFGS is to store the vector sequence in the calculation of approximation $H_{t+1}$, instead of storing complete matrix $H_{t}$. L-BFGS makes further consolidation for the update formula of $H_{t+1}$,

The above equation means that the inverse Hessian approximation $H_{t+1}$ can be obtained using the sequence pair $\{s_{l},u_{l}\}_{l=t-p+1}^{t}$. $H_{t+1}$ can be computed if we know pairs $\{s_{l},y_{l}\}_{l=t-p+1}^{t}$. In other words, instead of storing and calculating the complete matrix $H_{t+1}$, L-BFGS only computes the latest $p$ pairs of $\{s_{l},y_{l}\}$. According to the equation, a recursive procedure can be reached. When the latest $p$ steps are retained, the calculation of $H_{t+1}$ can be expressed as

The update direction $d_{t}=H_{t}g_{t}$ can be calculated, where $g_{t}$ is the gradient of the objective function $f$. The detailed algorithm is shown in Algorithms 3 and 4.

For more information about BFGS and L-BFGS algorithms, one can refer to . Recently, the batch L-BFGS on machine learning was proposed , which uses the overlapping mini-batches for consecutive samples for quasi-Newton update. It means that the calculation of $u_{t}$ becomes $u_{t}=\nabla_{S_{t+1}}f(\theta_{t+1})-\nabla_{S_{t}}f(\theta_{t}),$ where $S_{t}$ is a small subset of samples, meanwhile $S_{t+1}$ and $S_{t}$ are not independent, perhaps containing a relatively large overlap. Some numerical results in have shown that the modification in L-BFGS is effective in practice.

III-B3 Stochastic Quasi-Newton Method

In many large-scale machine learning models, it is necessary to use a stochastic approximation algorithm with each step of update based on a relatively small training subset . Stochastic algorithms often obtain the best generalization performances in large-scale learning systems . The quasi-Newton method only uses the first-order gradient information to approximate the Hessian matrix. It is a natural idea to combine the quasi-Newton method with the stochastic method, so that it can perform on large-scale problems. Online-BFGS and online-LBFGS are two variants of BFGS .

Consider the minimization of a convex stochastic function,

where $\xi$ is a random seed. We assume that $\xi$ represents a sample (or a set of samples) consisting of an input-output pair $(x,y)$. In machine learning $x$ typically represents an input and $y$ is the target output. $f$ usually has the following form:

where $h$ is a prediction model parameterized by $\theta$, and $l$ is a loss function. We define $f_{i}(\theta)=f(\theta;x_{i},y_{i})$, and use the empirical loss to define the objective,

Typically, if a large amount of training data is used to train the machine learning models, a better choice is using mini-batch stochastic gradient,

where subset $S_{t}\subset\{1,2,3\cdots N\}$ is randomly selected. $c$ is the cardinality of $S_{t}$ and $c\ll N$. Let $S_{t}^{H}\subset\{1,2,3,\cdots,N\}$ be a randomly chosen subset of the training samples and the stochastic Hessian estimate can be

where $c_{h}$ is the cardinality of $S_{t}^{H}$. With given stochastic gradient, a direct approach to develop stochastic quasi-Newton method is to transform deterministic gradients into stochastic gradients throughout the iterations, such as online-BFGS and online-LBFGS , which are two stochastic adaptations of the BFGS algorithms. Specifically, following the BFGS described in the previous section, $s_{t},u_{t}$ are modified as

One disadvantage of this method is that each iteration requires two gradient estimates. Besides this, a more worrying fact is that updating the inverse Hessian approximations in each step may not be reasonable . Then the stochastic quasi-Newton (SQN) method was proposed, which is to use sub-sampled Hessian-vector products to update $H_{t}$ by the LBFGS according to . Meanwhile, the authors proposed an effective approach that decouples the stochastic gradient and curvature estimate calculations to obtain a stable Hessian approximation. In particular, since

Based on these techniques, an SQN Framework was proposed, and the detailed procedure is shown in Algorithm 5.

In the above algorithm, $V=\{s_{t},u_{t}\}$ is a collection of $m$ displacement pairs, and $g_{t}$ is the current stochastic gradient $\nabla F_{S_{t}}(\theta_{t})$. Meanwhile, the matrix-vector product $H_{t}g_{t}$ can be computed by a two-loop recursion as described in the previous section. Recently, more and more work has achieved very good results in stochastic quasi-Newton. Specifically, a regularized stochastic BFGS method was proposed, which makes a corresponding analysis of the convergence of this optimization method . Further, an online L-BFGS was presented in . A linearly convergent method was proposed , which combines the L-BFGS method in with the variance reduction technique. Besides these, a variance reduced block L-BFGS method was proposed, which works by employing the actions of a sub-sampled Hessian on a set of random vectors .

To sum up, we have discussed the techniques of using stochastic methods in second-order optimization. The stochastic quasi-Newton method is a combination of the stochastic method and the quasi-Newton method, which makes the quasi-Newton method extend to large datasets. We have introduced the related work of the stochastic quasi-Newton method in recent years, which reflects the potential of the stochastic quasi-Newton method in machine learning applications.

III-B4 Hessian-Free Optimization Method

The main idea of Hessian-free (HF) method is similar to Newton’s method, which employs second-order gradient information. The difference is that the HF method is not necessary to directly calculate the Hessian matrix $H$. It estimates the product $Hv$ by some techniques, and thus is called “Hessian free”.

Consider a local quadratic approximation $Q_{\theta}(d_{t})$ of the object $F$ around parameter $\theta$,

where $d_{t}$ is the search direction. The HF method applies the conjugate gradient method to compute an approximate solution $d_{t}$ of the linear system,

where $B_{t}=H(\theta_{t})$ is the Hessian matrix, but in practice $B_{t}$ is often defined as $B_{t}=H(\theta_{t})+\lambda I,\ \lambda\geq 0$ . The new update is then given by

where $\eta_{t}$ is the step size that ensures sufficient decrease in the objective function, usually obtained by a linear search. According to , the basic framework of HF optimization is shown in Algorithm 6.

The advantage of using the conjugate gradient method is that it can calculate the Hessian-vector product without directly calculating the Hessian matrix. Because in the CG-algorithm, the Hessian matrix is paired with a vector, then we can compute the Hessian-vector product to avoid the calculation of the Hessian inverse matrix. There are many ways to calculate Hessian-vector products, one of which is calculated by a finite difference as

Sub-sampled Hessian-Free Method HF is a well-known method, and has been studied for decades in the optimization literatures, but has shortcomings when applied to deep neural networks with large-scale data . Therefore, a sub-sampled technique is employed in HF, resulting in an efficient HF method . The cost in each iteration can be reduced by using only a small sample set $S$ to calculate $Hv$. The objective function has the following form:

In the $t$th iteration, the stochastic gradient estimation can be written as

and the stochastic Hessian estimate is expressed as

As described above, we can obtain the approximate solution of direction $d_{t}$ by employing the CG method to solve the linear system,

in which the stochastic gradient and stochastic Hessian matrix are used. The basic framework of sub-sampled HF algorithm is given in .

A natural question is how to determine the size of $S^{H}_{t}$. On one hand, $S^{H}_{t}$ can be chosen small enough so that the total cost of CG iteration is not much greater than a gradient evaluation. On the other hand, $S^{H}_{t}$ should be large enough to get useful curvature information from Hessian-vector product. How to balance the size of $S^{H}_{t}$ is a challenge being studied .

III-B5 Natural Gradient

The natural gradient method can be potentially applied to any objective function which measures the performance of some statistical models . It enjoys richer theoretical properties when applied to objective functions based on the KL divergence between the model’s distribution and the target distribution, or certain approximation surrogates of these .

The traditional gradient descent algorithm is based on the Euclidean space. However, in many cases, the parameter space is not Euclidean, and it may have a Riemannian metric structure. In this case, the steepest direction of the objective function cannot be given by the ordinary gradient and should be given by the natural gradient .

We consider such a model distribution $p(y|x,\theta)$, and $\pi(x,y)$ is an empirical distribution. We need to fit the parameters $\theta\in R^{N}$. Assume that $x$ is an observation vector, and $y$ is its associated label. It has the objective function,

and we need to solve the optimization problem,

According to , the natural gradient can be transformed from a traditional gradient multiplied by a Fisher information matrix, i.e.,

where $F$ is the object function, $\bigtriangledown F$ is the traditional gradient, $\bigtriangledown_{N}F$ is the natural gradient, and $G$ is the Fisher information matrix, with the following form:

The update formula with the natural gradient is

We cannot ignore that the application of the natural gradient is very limited because of too much computation. It is expensive to estimate the Fisher information matrix and calculate its inverse matrix. To overcome this limitation, the truncated Newton’s method was developed , in which the inverse is calculated by an iterative procedure, thus avoiding the direct calculation of the inverse of the Fisher information matrix. In addition, the factorized natural gradient (FNG) and Kronecker-factored approximate curvature (K-FAC) methods were proposed to use the derivatives of probabilistic models to calculate the approximate natural gradient update.

III-B6 Trust Region Method

The update process of most methods introduced above can be described as $\theta_{t}+\eta_{t}d_{t}$. The displacement of the point in the direction of $d_{t}$ can be written as $s_{t}$. The typical trust region method (TRM) can be used for unconstrained nonlinear optimization problems , in which the displacement $s_{t}$ is directly determined without the search direction $d_{t}$.

For the problem $\min f_{\theta}(x)$, the TRM uses the second-order Taylor expansion to approximate the objective function $f_{\theta}(x)$, denoted as $q_{t}(s)$. Each search is done within the range of trust region with radius $\bigtriangleup_{t}$. This problem can be described as

where $g_{t}$ is the approximate gradient of the objective function $f(x)$ at the current iteration point $x_{t}$, $g_{t}\approx\nabla f(x_{t})$, $B_{t}$ is a symmetric matrix, which is the approximation of Hessian matrix $\nabla^{2}f_{\theta}(x_{t})$, and $\bigtriangleup_{t}>0$ is the radius of the trust region. If the $L_{2}$ norm is used in the constraint function, it becomes the Levenberg-Marquardt algorithm .

If $s_{t}$ is the solution of the trust region subproblem (80), the displacement $s_{t}$ of each update is limited by the trust region radius $\bigtriangleup_{t}$. The core part of the TRM is the update of $\bigtriangleup_{t}$. In each update process, the similarity of the quadratic model $q(s_{t})$ and the objective function $f_{\theta}(x)$ is measured, and $\bigtriangleup_{t}$ is updated dynamically. The actual amount of descent in the $t$th iteration is

The predicted drop in the $t$th iteration is

The ratio $r_{t}$ is defined to measure the approximate degree of both,

It indicates that the model is more realistic than expected when $r_{t}$ is close to 1, and then we should consider expanding $\bigtriangleup_{t}$. At the same time, it indicates that the model predicts a large drop and the actual drop is small when $r_{t}$ is close to 0, and then we should reduce $\bigtriangleup_{t}$. Moreover, if $r_{t}$ is between 0 and 1, we can leave $\bigtriangleup_{t}$ unchanged. The thresholds 0 and 1 are generally set as the left and right boundaries of $r_{t}$ .

III-B7 Summary

We summarize the mentioned high-order optimization methods in terms of properties, advantages and disadvantages in Table II.

III-C Derivative-Free Optimization

For some optimization problems in practical applications, the derivative of the objective function may not exist or is not easy to calculate. The solution of finding the optimal point, in this case, is called derivative-free optimization, which is a discipline of mathematical optimization . It can find the optimal solution without the gradient information.

There are mainly two types of ideas for derivative-free optimization. One is to use heuristic algorithms. It is characterized by empirical rules and chooses methods that have already worked well, rather than derives solutions systematically. There are many types of heuristic optimization methods, including classical simulated annealing arithmetic, genetic algorithms, ant colony algorithms, and particle swarm optimization . These heuristic methods usually yield approximate global optimal values, and theoretical support is weak. We do not focus on such techniques in this section. The other is to fit an appropriate function according to the samples of the objective function. This type of method usually attaches some constraints to the search space to derive the samples. Coordinate descent method is a typical derivative-free algorithm , and it can be extended and applied to optimization algorithms for machine learning problems easily. In this section, we mainly introduce the coordinate descent method.

The coordinate descent method is a derivative-free optimization algorithm for multi-variable functions. Its idea is that a one-dimensional search can be performed sequentially along each axis direction to obtain updated values for each dimension. This method is suitable for some problems in which the loss function is non-differentiable.

The vanilla approach is to select a set of bases $e_{1},e_{2},...,e_{D}$ in the linear space as the search directions and minimizes the value of the objective function in each direction. For the target function $L(\Theta)$, when $\Theta^{t}$ is already obtained, the $j$th dimension of $\Theta^{t+1}$ is solved by

Thus, $L(\Theta^{t+1})\leq L(\Theta^{t})\leq...\leq L(\Theta^{0})$ is guaranteed. The convergence of this method is similar to the gradient descent method. The order of update can be an arbitrary arrangement from $e_{1}$ to $e_{D}$ in each iteration. The descent direction can be generalized from the coordinate axis to the coordinate block .

The main difference between the coordinate descent and the gradient descent is that each update direction in the gradient descent method is determined by the gradient of the current position, which may not be parallel to any coordinate axis. In the coordinate descent method, the optimization direction is fixed from beginning to end. It does not need to calculate the gradient of the objective function. In each iteration, the update is only executed along the direction of one axis, and thus the calculation of the coordinate descent method is simple even for some complicated problems. For indivisible functions, the algorithm may not be able to find the optimal solution in a small number of iteration steps. An appropriate coordinate system can be used to accelerate the convergence. For example, the adaptive coordinate descent method takes principal component analysis to obtain a new coordinate system with as little correlation as possible between the coordinates . The coordinate descent method still has limitations when performing on the non-smooth objective function, which may fall into a non-stationary point.

III-D Preconditioning in Optimization

Preconditioning is a very important technique in optimization methods. Reasonable preconditioning can reduce the iteration number of optimization algorithms. For many important iterative methods, the convergence depends largely on the spectral properties of the coefficient matrix . It can be simply considered that the pretreatment is to transform a difficult linear system $A\theta=b$ into an equivalent system with the same solution but better spectral characteristics. For example, if $M$ is a nonsingular approximation of the coefficient matrix $A$, the transformed system,

will have the same solution as the system $A\theta=b$. But (85) may be easier to solve and the spectral properties of the coefficient matrix $M^{-1}A$ may be more favorable.

In most linear systems, e.g., $A\theta=b$, the matrix $A$ is often complex and makes it hard to solve the system. Therefore, some transformation is needed to simplify this system. $M$ is called the preconditioner. If the matrix after using preconditioner is obviously structured, or sparse, it will be beneficial to the calculation .

The conjugate gradient algorithm mentioned previously is the most commonly used optimization method with preconditioning technology, which speeds up the convergence. The algorithm is shown in Algorithm 7.

III-E Public Toolkits for Optimization

Fundamental optimization methods are applied in machine learning problems extensively. There are many integrated powerful toolkits. We summarize the existing common optimization toolkits and present them in Table LABEL:3.5.

IV Developments and Applications for Selected Machine Learning Fields

Optimization is one of the cores of machine learning. Many optimization methods are further developed in the face of different machine learning problems and specific application environments. The machine learning fields selected in this section mainly include deep neural networks, reinforcement learning, variational inference and Markov chain Monte Carlo.

The deep neural network (DNN) is a hot topic in the machine learning community in recent years. There are many optimization methods for DNNs. In this section, we introduce them from two aspects, i.e., first-order optimization methods and high-order optimization methods.

The stochastic gradient optimization method and its adaptive variants have been widely used in DNNs and have achieved good performance. SGD introduces the learning rate decay factor and AdaGrad accumulates all previous gradients so that their learning rates are continuously decreased and converge to zero. However, the learning rates of these two methods make the update slow in the later stage of optimization. AdaDelta, RMSProp, Adam and other methods use the exponential averaging to provide effective updates and simplify the calculation. These methods use exponential moving average to alleviate the problems caused by the rapid decay of the learning rate but limit the current learning rate to only relying on a few gradients . Reddi et al. used a simple convex optimization example to demonstrate that the RMSProp and Adam algorithms could not converge . Almost all the algorithms that rely on a fixed-size window of the past gradients will suffer from this problem, including AdaDelta and Nesterov-accelerated adaptive moment estimation (Nadam) .

It is better to rely on the long-term memory of past gradients rather than the exponential moving average of gradients to ensure convergence. A new version of Adam , called AmsGrad, uses a simple correction method to ensure the convergence of the model while preserving the original computational performance and advantages. Compared with the Adam method, the AmsGrad makes the following changes to the first-order moment estimation and the second-order moment estimation:

where $\beta_{1t}$ is a non-constant which decreases with time, and $\beta_{2}$ is a constant learning rate. The correction is operated in the second-order moment $V_{t}$, making $\hat{V}_{t}$ monotonous. $\hat{V}_{t}$ is substantially used in the iteration of the target function. The AmsGrad method takes the long-term memory of past gradients based on the Adam method, guarantees the convergence in the later stage, and works well in applications.

Further, adjusting parameters $\beta_{1},\beta_{2}$ at the same time helps to converge to a certain extent. For example, $\beta_{1}$ can decay modestly as $\beta_{1t}=\frac{\beta_{1}}{t}$, $\beta_{1t}\leq\beta_{1}$, for all $t\in[T]$. $\beta_{2}$ can be set as $\beta_{2t}=1-\frac{1}{t}$, for all $t\in[T]$, as in AdamNC algorithm .

Another idea of combining SGD and Adam was proposed for solving the non-convergence problem of adaptive gradient algorithm . Adaptive algorithms, such as Adam, converge fast and are suitable for processing sparse data. SGD with momentum can converge to more accurate results. The combination of SGD and Adam develops the advantages of both methods. Specifically, it first trains with Adam to quickly drop and then switches to SGD for precise optimization based on the previous parameters at an appropriate switch point. The strategy is named as switching from Adam to SGD (SWATS) . There are two core problems in SWATS. One is when to switch from Adam to SGD, the other is how to adjust the learning rate after switching the optimization algorithm. The SWATS approach is described in detail below.

The movement $d^{Adam}$ of the parameter at iteration $t$ of the Adam is

where $\eta^{Adam}$ is the learning rate of Adam . The movement $d^{SGD}$ of the parameter at iteration $t$ of the SGD is

where $\eta^{SGD}$ is the learning rate of SGD and $g_{t}$ is the gradient of the current position .

The movement of SGD can be decomposed into the learning rates along Adam’s direction and its orthogonal direction. If SGD is going to finish the trajectory but Adam has not finished due to the momentum after selecting the optimization direction, walking along Adam’s direction is a good choice for SWATS. At the same time, SWATS also adjusts its optimized trajectory by moving in the orthogonal direction. Let

where $Proj_{Adam}$ means the projection in the direction of Adam. To reduce noise, a moving average can be used to correct the estimate of the learning rate,

Recently some researchers are trying to explain and improve the adaptive methods . Their strategies can also be combined with the above switching techniques to enhance the performance of the algorithm.

General fully connected neural networks cannot process sequential data such as text and audio. Recurrent neural network (RNN) is a kind of neural networks that is more suitable for processing sequential data. It was generally considered that the use of first-order methods to optimize RNN was not effective, because the SGD and its variant methods were difficult to learn long-term dependencies in sequence problems .

In recent years, a well-designed method of random parameter initialization scheme using only SGD with momentum without curvature information has achieved good results in training RNNs . In , some techniques for improving the optimization in training RNNs are summarized such as the momentum methods and NAG. The first-order optimization methods have got development for training RNNs, but they still face the problem of slow convergence in deep RNNs. The high-order optimization methods employing curvature information can accelerate the convergence near the optimal value and is considered to be more effective in optimizing DNNs.

IV-A2 High-Order Gradient Method in Deep Neural Networks

We have described the first-order optimization method applied in DNNs. As most DNNs use large-scale data, different versions of stochastic gradient methods were developed and have got excellent performance and properties. For making full use of gradient information, the second-order method is gradually applied to DNNs. In this section, we mainly introduce the Hessian-free method in DNN.

Hessian-free (HF) method has been studied for a long time in the field of optimization, but it is not directly suitable for dealing with neural networks . As the objective function in DNN is not convex, the exact Hessian matrix may not be positive definite. Therefore, some modifications need to be made so that the HF method can be applied to neural networks .

The Generalized Gauss-Newton Matrix One solution is to use the generalized Gauss-Newton (GGN) matrix, which can be seen as an approximation of a Hessian matrix . The GGN matrix is a provably positive semidefinite matrix, which avoids the trouble of negative curvature. There are at least two ways to derive the GGN matrix . Both of them require that $f(\theta)$ can be expressed as a composition of two functions written as $f(\theta)=Q(F(\theta))$ where $f(\theta)$ is the object function and $Q$ is convex. The GGN matrix $G$ takes the following form,

Damping Methods Another modification to the HF method is to use different damping methods. For example, Tikhonov damping, one of the most famous damping methods, is implemented by introducing a quadratic penalty term into the quadratic model. A quadratic penalty term $\frac{\lambda}{2}d^{\top}d$ is added to the quadratic model,

where $B=H+\lambda I$, and $\lambda>0$ determines the “strength” of the damping which is a scalar parameter. Thus, $Bv$ is formulated as $Bv=(H+\lambda I)v=Hv+\lambda v$. However, the basic Tikhonov damping method is not good in training RNNs . Due to the complex structure of RNNs, the local quadratic approximation in certain directions in the parameter space, even at very small distances, maybe highly imprecise. The Tikhonov damping method can only compensate for this by increasing punishment in all directions because the method lacks a selective mechanism . Therefore, the structural damping was proposed, which makes the performance substantially better and more robust.

The HF method with structural damping can effectively train RNNs . Now we briefly introduce the HF method with structural damping. Let $e(x,\theta)$ mean the vector-value function of $\theta$ which can be interpreted as intermediate quantities during the calculation of $f(x,\theta)$, where $f(x,\theta)$ is the object function. For instance, $e(x,\theta)$ might contain the activation function of some layers of hidden units in neural networks (like RNNs). A structural damping can be defined as

where $D$ is a distance function or a loss function. It can prevent a large change in $e(x,\theta)$ by penalizing the distance between $e(x,\theta)$ and $e(x,\theta_{t})$. Then, the damped local objective can be written as

where $\mu$ and $\lambda$ are two parameters to be dynamically adjusted. $d$ is the direction at the $t$th iteration. More details of the structural damping can refer to .

Besides, there are many second-order optimization methods employed in RNNs. For example, quasi-Newton based optimization and L-BFGS were proposed to train RNNs .

In order to make the damping method based on punishment work better, the damping parameters can be adjusted continuously. A Levenberg-Marquardt style heuristic method was used to adjust $\lambda$ directly . The Levenberg-Marquardt heuristic is described as follows:

If $\gamma<\frac{1}{4}\lambda$ then $\lambda\leftarrow\frac{3}{2}\lambda$,

If $\gamma>\frac{3}{4}\lambda$ then $\lambda\leftarrow\frac{2}{3}\lambda$,

where $\gamma$ is a “reduction rate” with the following form,

Sub-sampling As sub-sampling Hessian can be used to handle large-scale data, several variations of the sub-sampling methods were proposed , which used either stochastic gradients or exact gradients. These approaches use $B_{t}=\nabla^{2}_{S_{t}}f(\theta_{t})$ as a Hessian approximation, where $S_{t}$ is a subset of samples. We need to compute the Hessian-vector product in some optimization methods. If we adopt the sub-sampling method, it also means that we can save a lot of computation in each iteration, such as the method proposed in .

Preconditioning Preconditioning can be used to simplify the optimization problems. For example, preconditioning can accelerate the CG method. It is found that diagonal matrices are particularly effective and one can use the following preconditioner :

where $\odot$ denotes the element-wise product and the exponent $\alpha$ is chosen to be less than 1.

IV-B Optimization in Reinforcement Learning

Reinforcement learning (RL) is an important research field of machine learning and is also one of the most popular topics. Agents using deep reinforcement learning have achieved great success in learning complex behavior skills and solved challenging control tasks in high-dimensional primitive perceptual state space . It interacts with the environment through the trial-and-error mechanism and learns optimal strategies by maximizing cumulative rewards .

We describe several concepts of reinforcement learning as follows:

Agent: making different actions according to the state of the external environment, and adjusting the strategy according to the reward of the external environment.

Environment: all things outside the agent that will be affected by the action of the agent. It can change the state and provide the reward to the agent.

State $s$: a description of the environment.

Action $a$: a description of the behavior of the agent.

Reward $r_{t}(s_{t-1},a_{t-1},s_{t})$: the timely return value at time $t$.

Policy $\pi(a|s)$: a function that the agent decides the action $a$ according to the current state $s$.

State transition probability $p(s^{\prime}|s,a)$: the probability distribution that the environment will transfer to state $s^{\prime}$ at the next moment, after the agent selecting an action $a$ based on the current state $s$.

$p(s^{\prime},r|s,a)$: the probability that the agent transforms to state $s^{\prime}$ and obtains the reward $r$, where the agent is in state $s$ and selecting the action $a$.

Many reinforcement learning problems can be described by Markov decision process (MDP) $<S,A,P,\gamma,r>$ , in which $S$ is state space, $A$ is action space, $P$ is state transition probability function, $r$ is reward function and $\gamma$ is the discount factor $0<\gamma<1$. At each time, the agent accepts a state and selects the action from an action set according to the policy. The agent receives feedback from the environment and then moves to the next state. The goal of reinforcement learning is to find a strategy that allows us to get the maximum $\gamma$-discounted cumulative reward. The discounted return is calculated by

People do not necessarily know the MDP behind the problem. From this point, reinforcement learning is divided into two categories. One is model-based reinforcement learning which knows the MDP of the whole model (including the transition probability $P$ and reward function $r$), and the other is the model-free method in which the MDP is unknown. Systematic exploration is required in the latter methods.

The most commonly used value function is the state value function,

which is the expected return of executing policy $\pi$ from state $s$. The state-action value function is also essential which is the expected return for selecting action $a$ under state $s$ and policy $\pi$,

The value function of the current state $s$ can be calculated by the value function of the next state $s^{\prime}$. The Bellman equations of $V_{\pi}(s)$ and $Q_{\pi}(s,a)$ describe the relation by

There are many reinforcement learning methods based on value function. They are called value-based methods, which play a significant role in RL. For example, Q-learning and SARSA are two popular methods which use temporal difference algorithms. The policy-based approach is to optimize the policy $\pi_{\theta}(a|s)$ directly and update the parameters $\theta$ by gradient descent .

The actor-critic algorithm is a reinforcement learning method combining policy gradient and temporal differential learning, which learns both a policy and a state value function. It estimates the parameters of two structures simultaneously.

The actor is a policy function, which is to learn a policy $\pi_{\theta}(a|s)$ to obtain the highest possible return.

The critic refers to the learned value function $V_{\phi}(s)$, which estimates the value function of the current policy, that is to evaluate the quality of the actor.

In the actor-critic method, the critic solves a problem of prediction, while the actor pays attention to the control . There is more information of actor-critic method in

The summary of the value-based method, the policy-based method, and the actor-critic method are as follows:

The value-based method: It needs to calculate value function, and usually gets a definite policy.

The policy-based method: It optimizes the policy $\pi$ without selecting an action according to value function.

The actor-critic method: It combines the above two methods, and learns both the policy $\pi$ and the state value function.

Deep reinforcement learning (DRL) combines reinforcement learning and deep learning, which defines problems and optimizes goals in the framework of RL, and solves problems such as state representation and strategy representation using deep learning techniques.

DRL has achieved great success in many challenging control tasks and uses DNNs to represent the control policy. For neural network training, a simple stochastic gradient algorithm or other first-order algorithms are usually chosen, but these algorithms are not efficient in exploring the weight space, which makes DRL methods often take several days to train . So, a distributed method was proposed to solve this problem, in which parallel actor-learners have a stabilizing effect during training . It executes multiple agents to interact with the environment simultaneously, which reduces the training time. But this method ignores the sampling efficiency. A scalable and sample-efficient natural gradient algorithm was proposed, which uses a Kronecker-factored approximation method to compute the natural policy gradient update, and employ the update to the actor and the critic (ACKTR) .

IV-C Optimization in Meta Learning

Meta learning is a popular research direction in the field of machine learning. It solves the problem of learning to learn. In the past cognition, the research of machine learning is to obtain a large amount of data in a specific task firstly and then use the data to train the model. In machine learning, adequate training data is the guarantee of achieving good performance. However, human beings can well process new tasks with only a few training samples, which are much more efficient than traditional machine learning methods. The key point could be that the human brain has learned “how to learn” and can make full use of past knowledge and experience to guide the learning of new tasks. Therefore, how to make machines have the ability to learn efficiently like human beings has become a frontier issue in machine learning.

The goal of meta learning is to design a model that can training well in the new tasks using as few samples as possible without overfitting. The process of adapting to the new tasks is essentially a learning process in the meta-testing, but only with limited samples from new tasks. The application of meta learning methods in supervised learning can solve the few-shot learning problems .

As few-shot learning problems receive more and more attention, meta learning is also developing rapidly. In general, meta learning methods can be summarized into the following three types : metric-based methods , model-based methods and optimization-based methods . In this subsection, we focus on the optimization-based meta learning methods. In meta learning, there are usually some tasks with sufficient training samples and a new task with only a few training samples. The main idea can be described as follows: in the meta-train step, sample a task $\tau$ from the total task set $\mathcal{T}$, which contains $(D_{\tau}^{train},D_{\tau}^{test})$. For task $\tau$, train and update the optimizer parameter $\theta$ with the training samples $D_{\tau}^{train}$, update the meta-optimizer parameter $\phi$ with the test samples $D_{\tau}^{test}$. The process of sampling tasks and updating parameters are repeated multiple times. In the meta-test step, the trained meta-optimizer is used for learning a new task.

Since the purpose of meta learning is to achieve fast learning, a key point is to make the gradient descent more accurately in the optimization. In some meta learning methods, the optimization process itself can be regarded as a learning problem to learn the prediction gradient rather than a determined gradient descent algorithm . Neural networks with original gradient as input and prediction gradient as output is often used as a meta-optimizer . The neural work is trained using the training and test samples from other tasks and used in the new task. The parameter update in the process of training is as follows:

where $\theta_{t}$ is the model parameter at the iteration $t$, and $N$ is the meta-optimizer with parameter $\phi$ that learns how to predict the gradient. After training, the meta-optimizer $N$ and its parameter $\phi$ are updated according to the loss value in the test samples. The experiments have confirmed that learning neural optimizers is advantageous compared to the most advanced adaptive stochastic gradient optimization methods used in deep learning . Due to the similarity between the gradient update in backpropagation and the cell state update in the long short-term memory (LSTM), LSTM is often used as the meta-optimizer .

A model-agnostic meta learning algorithm (MAML) is another method for meta learning which was proposed to learn the parameters of any model subjected to gradient descent methods. It is applicable to different learning problems, including classification, regression and reinforcement learning . The basic idea of the model-agnostic algorithm is to begin multiple tasks at the same time, and then get the synthetic gradient direction of different tasks, so as to learn a common base model. The main process can be described as follows: in the meta-train step, multiple tasks batch $\tau_{i}$, which contains $(D^{train}_{i},D^{test}_{i})$, are extracted from the total task set $\mathcal{T}$. For all $\tau_{i}$, train and update the parameter $\theta_{i}^{{}^{\prime}}$ with the train samples $D^{train}_{i}$:

where $\alpha$ is the learning rate of training process and $J_{\tau_{i}}(\theta)$ is the loss function in task $i$ with training samples $D_{i}^{train}$. After the training step, use the synthetic gradient direction of these parameters $\theta_{i}^{{}^{\prime}}$ on the test samples $D^{test}_{i}$ of the respective task to update parameter $\theta$:

where $\beta$ is the meta learning rate of the test process and $J_{\tau_{i}}(\theta)$ is the loss function in task $i$ with test samples $D_{i}^{test}$. The meta-train step is repeated multiple times to optimize a good initial parameter $\theta$. In the meta-test step, the trained parameter $\theta$ is used as the initial parameter such that the model has a maximal performance on the new task. MAML does not introduce additional parameters for meta learning, nor does it require a specific learner architecture. The development of the method is of great significance to the optimization-based meta learning methods. Recently, an expanded task-agnostic meta learning algorithm is proposed to enhance the generalization of meta-learner towards a variety of tasks, which achieves outstanding performance on few-shot classification and reinforcement learning tasks .

IV-D Optimization in Variational Inference

In the machine learning community, there are many attractive probabilistic models but with complex structures and intractable posteriors, and thus some approximate methods are used, such as variational inference and Markov chain Monte Carlo (MCMC) sampling. Variational inference, a common technique in machine learning, is widely used to approximate the posterior density of the Bayesian model, which transforms intricate inference problems into high-dimensional optimization problems . Compared with MCMC, the variational inference is faster and more suitable for dealing with large-scale data. Variational inference has been applied to large-scale machine learning tasks, such as large-scale document analysis, computer vision and computational neuroscience .

Variational inference often defines a flexible family of distributions indexed by free parameters on latent variables , and then finds the variational parameters by solving an optimization problem.

Now let us review the principle of variational inference . Variational inference approximates the true posterior by attempting to minimize the Kullback-Leibler (KL) divergence between a potential factorized distribution and the true posterior.

Let $Z=\{z_{i}\}$ represent the set of all latent variables and parameters in the model and $X=\{x_{i}\}$ be a set of all observed data. The joint likelihood of $X$ and $Z$ is $p(Z,X)=p(Z)p(X|Z)$. In Bayesian models, the posterior distribution $p(Z|X)$ should be computed to make further inference.

What we need to do is to approximate $p(Z|X)$ with the distribution $q(Z)$ that belongs to a constrained family of distributions. The goal is to make the two distributions as similar as possible. Variational inference chooses KL divergence to measure the difference between the two distributions, that is to minimize the KL divergence of $q(Z)$ and $p(Z|X)$. Here is the formula for the KL divergence between $q$ and $p$:

Variational inference can be treated as an optimization problem with the goal of minimizing the evidence lower bound. A direct method is to solve this optimization problem using the coordinate ascent, which is called coordinate ascent variational inference (CAVI). CAVI iteratively optimizes each factor of the mean-field variational density, while holding the others fixed .

Then the CAVI algorithm can be given below in Algorithm 8.

In traditional coordinate ascension algorithms, the efficiency of processing large data is very low, because each iteration needs to compute all the data, which is very time-consuming. Modern machine learning models often need to analyze and process large-scale data, which is difficult and costly. Stochastic optimization enables machine learning to be extended on massive data . This reminds us of an attractive technique to handle large data sets: stochastic optimization . By introducing stochastic optimization into variational inference, the stochastic variational inference (SVI) was proposed , in which the exponential family is taken as a typical example.

Gaussian process (GP) is an important machine learning method based on statistical learning and Bayesian theory. It is suitable for complex regression problems such as high dimensions, small samples, and nonlinearities. GP has the advantages of strong generalization ability, flexible non-parametric inference, and strong interpretability. However, the complexity and storage requirements of accurate solution for GP are high, which hinders the development of GP under large-scale data. The stochastic variational inference method introduced in this section can popularize variational inference on large-scale datasets, but it can only be applied to probabilistic models with factorized structures. For GPs whose observations are correlated with each other, the stochastic variational inference can be adapted by introducing the global inducing variables as variational variables . Specifically, the observations are assumed to be conditionally independent given the inducing variables and the variational distribution for the inducing variables is assumed to have an explicit form. Thus, the resulting GP model can be factorized in a necessary manner, enabling the stochastic variational inference. This method can also be easily extended to models with non-Gaussian likelihood or latent variable models based on GPs.

IV-E Optimization in Markov Chain Monte Carlo

Markov chain Monte Carlo (MCMC) is a class of sampling algorithms to simulate complex distributions that are difficult to sample directly. It is a practical tool for Bayesian posterior inference. The traditional and common MCMC algorithms include Gibbs sampling, slice sampling, Hamiltonian Monte Carlo (HMC) , Reimann manifold variants , and so on. These sampling methods are limited by the computational cost and are difficult to extend to large-scale data.This section takes HMC as an example to introduce the optimization in MCMC. The bottleneck of the HMC is that the gradient calculation is costly on large data sets.

We first introduce the derivation of HMC. Consider the random variable $\theta$, which can be sampled from the posterior distribution,

where $D$ is the set of observations, and $U$ is the potential energy function with the following formula:

In HMC , an independent auxiliary momentum variable $r$ is introduced from Hamiltonian dynamic. The Hamiltonian function and the joint distribution of $\theta$ and $r$ are described by

where $M$ denotes the mass matrix, and $K(r)$ is the kinetic energy function. The process of HMC sampling is derived by simulating the Hamiltonian dynamic system,

Hamiltonian dynamic describes the continuous motion of a particle. Hamiltonian equations are numerically approximated by the discretized leapfrog integrator for practical simulating . The update equations are as follows :

where $B(\theta)=\frac{1}{2}\epsilon V(\theta)$ is the diffusion matrix .

Since the discretization of the dynamical system introduces noise, the Metropolis-Hastings (MH) correction step should be done after the leapfrog step. These MH steps require expensive calculations overall data in each iteration. Beyond that, there is an incorrect stationary distribution in the stochastic gradient variant of HMC. Thus, Hamiltonian dynamic was further modified, which minimizes the effect of the additional noise, achieves the invariant distribution and eliminates MH steps . Specifically, a friction term is added to the dynamical process of momentum update:

The introduced friction term is helpful for decreasing total energy $H(\theta,r)$ and weakening the effects of noise in the momentum update phase. The dynamical system is also the type of second-order Langevin dynamics with friction in physics, which can explore efficiently and counteract the effect of the noisy gradients and thus no MH correction is required. This second-order Langevin dynamic MCMC method, called SGHMC, is used to deal with sampling problems on large data sets .

Moreover, HMC is highly sensitive to hyper-parameters, such as the path length (step number) $L$ and the step size $\epsilon$. If the hyper-parameters are not set properly, the efficiency of the HMC will drop dramatically. There are some methods to optimize these two hyper-parameters instead of manually setting them.

The value of path length $L$ has a great influence on the performance of HMC. If $L$ is too small, the distance between the resulting sample points will be very close; if $L$ is too large, the resulting sample points will loop back, resulting in wasted computation. In general, manually setting $L$ cannot maximize the sampling efficiency of the HMC.

Matthew et al. proposed an extension of the HMC method called the No-U-Turn sampler (NUTS), which uses a recursive algorithm to generate a set of possible independent samples efficiently, and stops the simulation by discriminating the backtracking automatically. There is no need to set the step parameter $L$ manually. In models with multiple discrete variables, the ability of NUTS to select the track length automatically allows it to generate more valid samples and perform more efficiently than the original HMC.

IV-E2 Adaptive Step Size ϵitalic-ϵ\epsilon

The performance of HMC is highly sensitive to the step size $\epsilon$ in leapfrog integrator. If $\epsilon$ is too small, the update will slow, and the calculation cost will be high; if $\epsilon$ is too large, the rejection rate will be high, resulting in useless updates.

To set $\epsilon$ reasonably and adaptively, a vanishing adaptation of the dual averaging algorithm can be used in HMC . Specifically, a statistic $H_{t}=\delta-\alpha_{t}$ is adopted in dual averaging method, where $\delta$ is the desired average acceptance probability, and $\alpha_{t}$ is the current Metropolis-Hasting acceptance probability for iteration $t$. The statistic $H_{t}$’s expectation $h(\epsilon)$ is defined as

The hyper-parameters in the HMC include not only the step size $\epsilon$ and the length of iteration steps $L$, but also the mass $M$, etc. Optimizing these hyper-parameters can help improve sampling performance . It is convenient and efficient to tune the hyper-parameters automatically without cumbersome adjustments based on data and variables in MCMC. These adaptive tuning methods can be applied to other MCMC algorithms to improve the performance of the samplers.

In addition to second-order SGHMC, stochastic gradient Langevin dynamics (SGLD) is a first-order Langevin dynamic technique combined with stochastic optimization. Efficient variants of both SGLD and SGHMC are still active .

V Challenges and Open Problems

With the rise of practical demand and the increase of the complexity of machine learning models, the optimization methods in machine learning still face challenges. In this part, we discuss open problems and challenges for some optimization methods in machine learning, which may offer suggestions or ideas for future research and promote the wider application of optimization methods in machine learning.

There are still many challenges while optimizing DNNs. Here we mainly discuss two challenges with respect to data and model, respectively. One is insufficient data in training, and the other is a non-convex objective in DNNs.

In general, deep learning is based on big data sets and complex models. It requires a large number of training samples to achieve good training effects. But in some particular fields, finding a sufficient amount of training data is difficult. If we do not have enough data to estimate the parameters in the neural networks, it may lead to high variance and overfitting.

There are some techniques in neural networks that can be used to reduce the variance. Adding $L_{2}$ regularization to the objective is a natural method to reduce the model complexity. Recently, a common method is dropout . In the training process, each neuron is allowed to stop working with a probability of $p$, which can prevent the synergy between certain neurons. $M$ subnets can be sampled like bagging by multiple puts and returns . Each expected result at the output layer is calculated as

where $p(M_{i})$ is the probability of the $i$th subnet. Dropout can prevent overfitting and improve the generalization ability of the network, but its disadvantage is increasing the training time as each training changes from the full network to a sub-network .

Not only overfitting but also some training details will affect the performance of the model due to the complexity of the DNNs. The improper selection of the learning rate and the number of iterations in the SGD will make the model unable to converge, which makes the accuracy of model fluctuate greatly. Besides, taking an inappropriate black box of neural network construction may result in training not being able to continue, so designing an appropriate neural network model is particularly important. These impacts are even greater when data are insufficient.

The technology of transfer learning can be applied to build networks in the scenario of insufficient data. Its idea is that the models trained from other data sources can be reused in similar target fields after certain modifications and improvements, which dramatically alleviates the problems caused by insufficient datasets. Moreover, the advantages brought by transfer learning are not limited to reducing the need for sufficient training data, but also can avoid overfitting effectively and achieve better performance in general. However, if target data is not as relevant to the original training data, the transferred model does not bring good performance.

Meta learning methods can be used for systematically learning parameter initialization, which ensures that training begins with a suitable initial model. However, it is necessary to ensure the correlation between multiple tasks for meta-training and tasks for meta-testing. Under the premise of models with similar data sources for training, transfer learning and meta learning can overcome the difficulties caused by insufficient training data in new data sources, but these methods usually introduce a large number of parameters or complex parameter adjustment mechanisms, which need to be further improved for specific problems. Therefore, using insufficient data for training DNNs is still a challenge.

V-A2 Non-convex Optimization in Deep Neural Network

Convex optimization has good properties and a comprehensive set of tools are open to solve the optimization problem. However, many machine learning problems are formulated as non-convex optimization problems. For example, almost all the optimization problems in DNNs are non-convex. Non-convex optimization is one of the difficulties in the optimization problem. Unlike convex optimization, there may be innumerable optimum solutions in its feasible domain in non-convex problems. The complexity of the algorithm for searching the global optimal value is NP-hard .

In recent years, non-convex optimization has gradually attracted the attention of researches. The methods for solving non-convex optimization problems can be roughly divided into two types. One is to transform the non-convex optimization into a convex optimization problem, and then use the convex optimization method. The other is to use some special optimization method for solving non-convex functions directly. There is some work on summarizing the optimization methods for solving non-convex functions from the perspective of machine learning .

Relaxation method: Relax the problem to make it become a convex optimization problem. There are many relaxation techniques, for example, the branch-and-bound method called $\alpha$BB convex relaxation , which uses a convex relaxation at each step to compute the lower bound in the region. The convex relaxation method has been used in many fields. In the field of computer vision, a convex relaxation method was proposed to calculate minimal partitions . For unsupervised and semi-supervised learning, the convex relaxation method was used for solving semidefinite programming .

Non-convex optimization methods: These methods include projection gradient descent , alternating minimization , expectation maximization algorithm and stochastic optimization and its variants .

V-B Difficulties in Sequential Models with Large-Scale Data

When dealing with large-scale time series, the usual solutions are using stochastic optimization, processing data in mini-batches, or utilizing distributed computing to improve computational efficiency . For a sequential model, segmenting the sequences can affect the dependencies between the data on the adjacent time indices. If sequence length is not an integral multiple of the mini-batch size, the general operation is to add some items sampled from the previous data into the last subsequence. This operation will introduce the wrong dependency in the training model. Therefore, the analysis of the difference between the approximated solution obtained and the exact solution is a direction worth exploring.

Particularly, in RNNs, the problem of gradient vanishing and gradient explosion is also prone to occur. So far, it is generally solved by specific interaction modes of LSTM and GRU or gradient clipping. Better appropriate solutions for dealing with problems in RNNs are still worth investigating.

V-C High-Order Methods for Stochastic Variational Inference

The high-order optimization method utilizes curvature information and thus converges fast. Although computing and storing the Hessian matrices are difficult, with the development of research, the calculation of the Hessian matrix has made great progress , and the second-order optimization method has become more and more attractive. Recently, stochastic methods have also been introduced into the second-order method, which extends the second order method to large-scale data .

We have introduced some work on stochastic variational inference. It introduces the stochastic method into variational inference, which is an interesting and meaningful combination. This makes variational inference be able to handle large-scale data. A natural idea is whether we can incorporate second-order optimization methods (or higher-order) into stochastic variational inference, which is interesting and challenging.

V-D Stochastic Optimization in Conjugate Gradient

Stochastic methods exhibit powerful capabilities when dealing with large-scale data, especially for first-order optimization . Then the relevant experts and scholars also introduced this stochastic idea to the second-order optimization methods and achieved good results.

Conjugate gradient method is an elegant and attractive algorithm, which has the advantages of both the first-order and second-order optimization methods. The standard form of a conjugate gradient is not suitable for a stochastic approximation. Through using the fast Hessian-gradient product, the stochastic method is also introduced to conjugate gradient, in which some numerical results show the validity of the algorithm . Another version of stochastic conjugate gradient method employs the variance reduction technique, and converges quickly with just a few iterations and requires less storage space during the running process . The stochastic version of conjugate gradient is a potential optimization method and is still worth studying.

VI Conclusion

This paper introduces and summarizes the frequently used optimization methods from the perspective of machine learning, and studies their applications in various fields of machine learning. Firstly, we describe the theoretical basis of optimization methods from the first-order, high-order, and derivative-free aspects, as well as the research progress in recent years. Then we describe the applications of the optimization methods in different machine learning scenarios and the approaches to improve their performance. Finally, we discuss some challenges and open problems in machine learning optimization methods.

References