Stochastic Training of Neural Networks via Successive Convex Approximations

TL;DR

This paper introduces a novel stochastic successive convex approximation algorithm for neural network training that leverages first-order information, leading to faster convergence and better minima, with easy parallelization.

Contribution

The paper presents a new SCA-based training method for neural networks that uses only first-order stochastic approximations, improving convergence speed and scalability.

Findings

Outperforms state-of-the-art training techniques in convergence speed.

Effectively parallelizable across multiple computational units.

Demonstrates success on medium-sized and large-scale datasets.

Abstract

This paper proposes a new family of algorithms for training neural networks (NNs). These are based on recent developments in the field of non-convex optimization, going under the general name of successive convex approximation (SCA) techniques. The basic idea is to iteratively replace the original (non-convex, highly dimensional) learning problem with a sequence of (strongly convex) approximations, which are both accurate and simple to optimize. Differently from similar ideas (e.g., quasi-Newton algorithms), the approximations can be constructed using only first-order information of the neural network function, in a stochastic fashion, while exploiting the overall structure of the learning problem for a faster convergence. We discuss several use cases, based on different choices for the loss function (e.g., squared loss and cross-entropy loss), and for the regularization of the NN's weights. We experiment on several medium-sized benchmark problems, and on a large-scale dataset involving simulated physical data. The results show how the algorithm outperforms state-of-the-art techniques, providing faster convergence to a better minimum. Additionally, we show how the algorithm can be easily parallelized over multiple computational units without hindering its performance. In particular, each computational unit can optimize a tailored surrogate function defined on a randomly assigned subset of the input variables, whose dimension can be selected depending entirely on the available computational power.

Tables

Table 1. TABLE II: Schematic description of the mid-sized datasets. For the NN topology, x / y 𝑥 𝑦 x/y denotes a NN with two layers of dimensions x 𝑥 x and y 𝑦 y respectively.

Dataset	Samples	Features	NN Topology	Prediction task	Source
CASP	$45730$	$9$	$9/10/6/1$	Protein’s structure prediction from physicochemical properties	UCI
Parkinsons	$5875$	$19$	$19/15/5/1$	UPDRS score from measurements of a remote telemonitoring device	UCI
SkillCraft1	$3395$	$20$	$18/15/10/1$	Predicted level from in-game statistics	UCI
Wine	$4898$	$12$	$11/10/4/1$	Wine quality from chemical measurements	UCI

Table 2. TABLE III: Results (in terms of mean-squared error on the test set) of different algorithms. Best results for both groups are highlighted in bold.

Dataset	SGD	AdaGrad	RMSProp	Adam	Proposed
CASP	$0.2713\pm 0.0053$	$0.2348\pm 0.0074$	$0.2182\pm 0.0114$	$0.2107\pm 0.0103$	$0.2017\pm 0.0045$
Parkinsons	$0.1550\pm 0.0062$	$0.1389\pm 0.0042$	$0.1417\pm 0.0075$	$0.1403\pm 0.0058$	$0.1374\pm 0.0045$
SkillCraft1	$0.0708\pm 0.0040$	$0.0688\pm 0.0029$	$0.0738\pm 0.0083$	$0.0707\pm 0.0061$	$0.0675\pm 0.0034$
Wine	$0.0692\pm 0.0029$	$0.0551\pm 0.0024$	$0.0562\pm 0.0025$	$0.0543\pm 0.0034$	$0.0528\pm 0.0023$

Equations

$$\left\lVert\nabla f(\mathbf{w}_{1};\mathbf{x})-\nabla f(\mathbf{w}_{2};\mathbf{x})\right\rVert_{2}\leq L\bigl{\lVert}\mathbf{w}_{1}-\mathbf{w}_{2}\bigr{\rVert}_{2}\,,$$

$$\underset{\mathbf{w}}{\min}\;\;U(\mathbf{w})=\frac{1}{N}\sum_{i=1}^{N}l\bigl{(}y_{i},f(\mathbf{w};\mathbf{x}_{i})\bigr{)}+\lambda\cdot r(\mathbf{w})\,,$$

$${\mathbf{w}}_{n+1}={\mathbf{w}}_n-{\alpha }_n\tilde{\nabla }U({\mathbf{w}}_n),$$

$${\hat{\mathbf{w}}}_{n+1}=$$

$$\displaystyle\Biggl{.}+\left(1-\rho_{n}\right)\mathbf{d}_{n}^{T}\left(\mathbf{w}-\mathbf{w}_{n}\right)+\tau\left\lVert\mathbf{w}-\mathbf{w}_{n}\right\rVert_{2}^{2}\Biggr{\}}\,,$$

$${\mathbf{w}}_{n+1}=\left(1-{\alpha }_n\right){\mathbf{w}}_n+{\alpha }_n{\hat{\mathbf{w}}}_{n+1},$$

$$\mathbf{d}_{n+1}=\left(1-\rho_{n}\right)\mathbf{d}_{n}+\rho_{n}\cdot\biggl{(}\frac{1}{L}\sum_{i\in\mathcal{I}_{n}}\nabla l_{i}(\mathbf{w}_{n})\biggr{)}\,,$$

$$\underset{n\to \infty }{\lim }{\alpha }_n=0,\sum _{n=1}^{\infty }{\alpha }_n=\infty ,\sum _{n=1}^{\infty }{\alpha }_n^2<\infty .$$

$$\underset{n\to \infty }{\lim }{\rho }_n=0,\sum _{n=1}^{\infty }{\rho }_n=\infty ,\sum _{n=1}^{\infty }{\rho }_n^2<\infty .$$

$$\underset{n\to \infty }{\lim }{\alpha }_n/{\rho }_n=0.$$

$$\widetilde{l}_{i}(\mathbf{w};\mathbf{w}_{n})=l_{i}(\mathbf{w}_{n})+\nabla l_{i}(\mathbf{w}_{n})\bigl{(}\mathbf{w}-\mathbf{w}_{n}\bigr{)}\,.$$

$$\widehat{\mathbf{w}}_{n+1}=\frac{1}{\lambda}\biggl{(}\rho_{n}\frac{1}{L}\nabla\sum_{i\in\mathcal{I}_{n}}l_{i}(\mathbf{w}_{n})+\left(1-\rho_{n}\right)\mathbf{d}_{n}\biggr{)}\,.$$

$${\tilde{f}}_i(\mathbf{w};{\mathbf{w}}_n)=f({\mathbf{w}}_n;{\mathbf{x}}_i)+{\mathbf{J}}_{i,n}^T\left(\mathbf{w}-{\mathbf{w}}_n\right),$$

$$\widetilde{l}_{i}(\mathbf{w};\mathbf{w}_{n})=l\bigl{(}y_{i},\widetilde{f}_{i}(\mathbf{w};\mathbf{w}_{n})\bigr{)}\,.$$

$$l(a,b)\triangleq {\left(a-b\right)}^2,r(\mathbf{w})\triangleq \frac{1}{2}{\Vert \mathbf{w}\Vert }_2^2.$$

$$r_{i,n}=y_i-f({\mathbf{w}}_n;{\mathbf{x}}_i)+{\mathbf{J}}_{i,n}^T{\mathbf{w}}_n.$$

$$\widehat{\mathbf{w}}_{n+1}=\underset{\mathbf{w}}{\arg\min}\;\;\biggl{\{}\mathbf{w}^{T}\Bigl{(}\mathbf{A}_{n}+\lambda\mathbf{I}\Bigr{)}\mathbf{w}_{i}-2\mathbf{b}_{n}^{T}\mathbf{w}\biggr{\}}\,,$$

$${\mathbf{A}}_n$$

$${\mathbf{b}}_n$$

$$\widehat{\mathbf{w}}_{n+1}=\Bigl{(}\mathbf{A}_{n}+\lambda\mathbf{I}\Bigr{)}^{-1}\mathbf{b}_{n}\,.$$

$$\widehat{\mathbf{w}}_{n+1}=\Bigl{(}\mathbf{A}_{n}+\left(\lambda+\tau\right)\mathbf{I}\Bigr{)}^{-1}\left(\mathbf{b}_{n}+\tau\mathbf{w}_{n}\right)\,,$$

$$r(\mathbf{w})\triangleq {\Vert \mathbf{w}\Vert }_1=\sum _{i=1}^Q|w_i|.$$

$$\widehat{\mathbf{w}}_{n+1}=\underset{\mathbf{w}}{\arg\min}\;\;\biggl{\{}\mathbf{w}^{T}\mathbf{A}_{n}\mathbf{w}_{i}-2\mathbf{b}_{n}^{T}\mathbf{w}+\lambda\left\lVert\mathbf{w}\right\rVert_{1}\biggr{\}}\,.$$

$$\mathbf{w}=\bigcup _{p=1}^P{\mathbf{w}}_p.$$

$$r(\mathbf{w})\triangleq \sum _{p=1}^P{a}_p{\Vert {\mathbf{w}}_p\Vert }_2,$$

$$l(a,b)\triangleq a\log (b)+(1-a)\log (1-b).$$

$$f(\mathbf{w})=\sigma\bigl{(}f^{L}(\mathbf{w})\bigr{)}\,,$$

$$\widetilde{f}^{L}_{i}(\mathbf{w};\mathbf{w}_{n})=f^{L}(\mathbf{w}_{n};\mathbf{x}_{i})+\bigl{(}\mathbf{J}^{L}_{i,n}\bigr{)}^{T}\left(\mathbf{w}-\mathbf{w}_{n}\right)\,,$$

$$\widetilde{l}_{i}(\mathbf{w};\mathbf{w}_{n})=l_{i}\bigl{(}y_{i},\sigma(\widetilde{f}^{L}(\mathbf{w};\mathbf{w}_{n}))\bigr{)}\,.$$

$$r(\mathbf{w})=\sum _{i=1}^N\frac{1}{k}\sum _{j\in {\mathcal{N}}_i}{q}_{ij}{\Vert f(\mathbf{w};{\mathbf{x}}_i)-f(\mathbf{w};{\mathbf{x}}_j)\Vert }_2^2.$$

$$\tilde{r}(\mathbf{w};{\mathbf{w}}_n)=\sum _{i=1}^N\frac{1}{k}\sum _{j\in {\mathcal{N}}_i}{q}_{ij}{\Vert {\Delta }_{ij}-{\mathbf{J}}_{i,j,n}^T\mathbf{w}\Vert }_2^2,$$

$${\Delta }_{ij}$$

Full text

I Introduction

Training a neural network (NN) involves the minimization of a high-dimensional, non-convex loss function, whose gradients can further exponentially vanish or explode during training [1]. Until very recently, this combination appeared too challenging to be confronted directly, and the first wave of deep NNs were trained with expensive, layer-wise, greedy initialization procedures. Today, the general consensus is that poor minima are less of a problem than what previously believed (as opposed to, e.g., pathological saddle points [2]). Optimization tools have equivalently shifted to fully supervised routines, employing stochastic first-order algorithms, such as Adam [3], supplemented by advanced regularization methods, including dropout, batch normalization, and several others. As a result, stochastic gradient descent (SGD) has become a fundamental mainstay in machine learning [4].

To speed up the training process, several researchers have considered the inclusion of curvature information in the optimization step [5, 6, 7, 8], when working with mini-batches of the full training set. This type of second-order methods (e.g., quasi-Newton) looks for a descent direction by minimizing a quadratic approximation to the cost function [9]. In a stochastic setting, the relative Hessian information is estimated with a noisy version computed from the current batch of data points. However, this quadratic approximation can be unreliable, with the error due to stochasticity compounding and making the overall process hard to implement in a straightforward fashion. As a consequence, these algorithms have failed to gain wide recognition in subsequent years.

This dichotomy between SGD and second-order methods is not limited to the training of NNs, and several methods have been proposed in the optimization literature to overcome it. A very recent proposal in this sense is called the successive convex approximation (SCA) framework [10], which was initially proposed in the context of multi-agent systems [11]. Like second-order methods, SCA algorithms work by solving a series of approximations of the original optimization task, called surrogate problems. Differently from Newton techniques, however, these surrogates are not limited to simple quadratic approximations. In fact, any knowledge about the structure of the objective function can be exploited in order to design surrogates which as a result are (in general) both simple to compute and efficient to solve, as long as they maintain the first-order information of the original function (see Section III-B for more details). Recently, these techniques were also extended to a stochastic setting [12].

In light of the above, the aim of this paper is to explore stochastic SCA algorithms for training NN models. Specifically, we show that the structure of the optimization problem can be exploited to design particularly efficient surrogate problems, which are obtained by combining a linearization of the original NN model with the (convex) loss function, e.g., the squared loss or the cross-entropy loss. The resulting surrogates can be computed efficiently starting from mini-batches of data, without the need to compute expensive second-order information of the cost function. We focus specifically on the case of the squared loss, giving rise to quadratic optimization problems that can be solved immediately in closed form. As shown in the experimental section, they are able to provide excellent descent directions, surpassing most state-of-the-art stochastic solvers, including Adam [3] and Adagrad [13]. This is particularly useful in situations where we need fast convergence in a small number of iterations, such as in federated learning environments [14].

Contribution of the paper

We introduce a general framework for optimizing NN models using stochastic SCA algorithms, and we customize it for different loss functions, such as the squared loss, and different regularizers, including $\ell_{2}$ and $\ell_{1}$ norms to promote sparsity. We also consider the case of a non-convex regularization term, and show an immediate extension of the framework to handle it. In addition, we further build on the theory of SCA techniques [10] to provide a principled way to parallelize the computation (exploiting, e.g., a multi-core architecture), by defining surrogate problems on subsets of the overall vector of parameters, up to one weight per processor. This is useful for NN models comprising a large set of parameters, where different processors are in charge of optimizing separate portions of the network.

Outline of the paper

We analyze related literature in Section II. Next, Section III provides a general stochastic SCA algorithm and its convergence properties. In Section IV, we show how to customize the algorithm to the problem of NN training, and we explore several specializations employing different loss functions and/or regularizers. In Section V we describe how to parallelize the algorithm. Finally, we provide a comprehensive evaluation in Section VI, before concluding in Section VII.

Notation

We denote vectors using boldface lowercase letters, e.g., $\mathbf{a}$; matrices are denoted by boldface uppercase letters, e.g., $\mathbf{A}$. All vectors are assumed to be column vectors. The operator $\left\lVert\cdot\right\rVert_{p}$ is the standard $\ell_{p}$ norm on an Euclidean space. For $p=2$, it coincides with the Euclidean norm, while for $p=1$ we obtain the Manhattan (or taxicab) norm defined for a generic vector $\mathbf{v}\in\R^{B}$ as $\left\lVert\mathbf{v}\right\rVert_{1}=\sum_{k=1}^{B}|v_{k}|$. Other notation is introduced along the paper when required.

II Related works

The idea of successively replacing a non-convex objective with a series of convex approximations is not novel in the optimization literature, and it appears in a wide range of previous approaches, including convex-concave procedures [15] and proximal minimization algorithms [16]. However, most previous methods imposed stringent conditions on the approximant, such as it being a global upper bound of the original cost (e.g., the SUM algorithm in [17]). The SCA methods that we consider here originated in the context of multi-agent systems [11], and were later extended to deal with general non-convex optimization problems [18, 10, 19]. Under this framework, the (convex) approximation is only required to keep the first-order information of the original (non-convex) cost with respect to the current estimate, thus making its definition highly flexible. Additionally, the optimization problems can be easily decomposed into subproblems (see later in Section V), and convergence to a stationary point can be guaranteed under mild conditions. Several extensions were made to the basic framework, most notably SCA techniques for decentralized environments [20, 21], asynchronous processors [22], and stochastic updates [12].

To the best of our knowledge, the only works that applied SCA techniques for training NNs are [23, 24]. However, both are specific to a distributed setting with full batch updates, which is not scalable to the training of NNs with a large-scale dataset or with many parameters. The present paper significantly extends [24] to the case of stochastic updates computed from mini-batches of the training data. More tangentially related to our paper are the investigations in [25] and [26], which applied SCA techniques for training support vector models, always in a batch setting.

More generally, our definition of the surrogate function (see Section IV-A) builds on the fact that the non-convexity is specific to the NN model, while the loss and regularization terms are generally convex with respect to the NN output. Historically, the first work to exploit this idea is [27]. In [27], it is shown that training a NN with a single-hidden layer, where the number of neurons is also learned, is a convex problem with respect to the weights. Using a sparsity penalty, it is possible to replace the original training problem with a series of convex problems where a single hidden unit is incrementally added. Similar ideas are explored in [28]. Another paper exploiting a related idea is the ‘mollifying network’ presented in [29], where the original problem is successively relaxed by convolving it with a mollifier function. Convexity is guaranteed only in the first iteration of training. None of these papers, however, is connected to the SCA techniques described next.

III Preliminaries

We begin by stating the NN optimization problem in Section III-A. Then, we describe a generic stochastic SCA algorithm for non-convex optimization problems in Section III-B. As stated before, in order to be applied SCA techniques require the definition of a proper surrogate cost function. How to exploit the structure of our optimization problem to this end, together with several practical use cases, is the topic of the next section.

III-A Formulation of the problem

We want to train a generic NN $f(\mathbf{w};\mathbf{x})$, which takes as input a real-valued vector $\mathbf{x}\in\R^{d}$, and outputs a scalar value. We consider networks with a single output mostly for readability, but everything can be extended seamlessly to the case of multiple outputs. The output of the network depends on a set of $Q$ adaptable parameters (e.g., the weights connecting the layers), that we collect in a single vector $\mathbf{w}\in\R^{Q}$ to be optimized depending on some training data. Additionally, note that we consider the NN as a function of its parameters, in order to make the following notation simpler.

The NN can have any number of layers, nonlinearities, etc., as long as the following assumption is satisfied for any possible input vector.

Assumption A [On the NN model]:

(A1)

$f\in C^{1}$, i.e., it is continuously differentiable with respect to $\mathbf{w}$;

(A2)

$f$ has Lipschitz continuous gradient, with respect to $\mathbf{w}$, for some Lipschitz constant $L>0$, i.e.:111Unless specified, all gradients in the manuscript are taken with respect to $\mathbf{w}$ or its current estimate, which is always clear from the context.

[TABLE]

for any $\mathbf{w}_{1},\mathbf{w}_{2}\in\R^{Q}$.

These are satisfied by most NN models currently used in the literature, the most notable exception being activation functions with a finite number of non-differentiable points such as ReLu neurons [30], maxout neurons [31], piecewise linear adaptable functions [32], and a few others. In fact, the aforementioned cases lead to non-convex, non-differentiable neural networks functions that are not currently handled by the SCA method. An interesting future line of research will be to design optimization methods specifically tailored for such important cases.

We are provided with a training dataset of $N$ input/output pairs $\mathcal{S}=\left\{\mathbf{x}_{i},y_{i}\right\}$, and the learning task aims at solving the following regularized optimization problem:

[TABLE]

where $l(\cdot,\cdot)$ is a convex, smooth loss function also satisfying the conditions in Assumption A, while $r(\cdot)$ is a (possibly non-smooth) convex regularization term, and $\lambda>0$ is a user-defined scalar that weights the two terms. Typical losses used in training NNs are the squared error for regression and the cross-entropy loss for classification, while typical regularizers are $\ell_{1}/\ell_{2}$ penalties or a combination of them. All loss functions used in practice are convex with respect to their arguments, with the non-convexity of $f$ making the overall problem non-convex. This point will be essential for the development of the algorithm in the following.

A classical approach to solve (2) is to use a stochastic, first-order optimization algorithm, where at the $n$th iteration we sample $L$ indexes in $\left\{1,\ldots,N\right\}$, with $L\ll N$. We denote by $\mathcal{I}_{n}$ the random variable containing the indexes sampled at time $n$, and by $\mathcal{B}_{n}$ the corresponding mini-batch of elements extracted from our dataset. We update the current estimate $\mathbf{w}_{n}$ following (noisy) gradient information from $U(\mathbf{w})$ [13]. In particular, using $\mathcal{B}_{n}$ to compute a noisy gradient $\widetilde{\nabla}U(\mathbf{w}_{n})$ of (2), the simplest possibility is stochastic gradient descent (SGD):

[TABLE]

where $\alpha_{n}$ is a (generally decreasing) sequence of step-sizes. This basic strategy can be accelerated in a number of ways, including weight-dependent step-sizes, momentum, gradient averaging, and so on. We refer to [13] for an up-to-date survey on the topic. See [33] for a recent critique on some of these accelerated methods for training non-convex models.

III-B Stochastic SCA optimization

By only exploiting first-order information on the cost function $U(\mathbf{w})$, SGD and its variants can incur in slow convergence speed and, more in general, they do not leverage efficiently all the information contained in the mini-batch $\mathcal{B}_{n}$ [6]. In contrast, Newton and quasi-Newton methods try to add curvature information to the optimization process, by iteratively minimizing a (noisy) quadratic approximation of the cost function. However, as we discussed in the introduction, these methods require several adjustments to be efficient while training NNs [6, 34], and they have failed to gain widespread adoption. SCA techniques try to surpass these disadvantages by building approximations of the cost function that aim to preserve as much as possible of its ‘hidden convexity’. In this section we describe the general stochastic SCA technique introduced in [12], which is used as the building block for the optimization algorithms of the next section.

Simplifying the notation, denote by l_{i}(\mathbf{w})=l\bigl{(}y_{i},f(\mathbf{w};\mathbf{x}_{i})\bigr{)} the $i$th loss term in (2) as a function of the NN parameters. Let us consider a surrogate loss $\widetilde{l}_{i}(\mathbf{w};\mathbf{w}_{n})$ of $l_{i}$ (with respect to the current weights’ estimate), having the following properties:

Assumption B [On the surrogate function]:

(B1)

$\widetilde{l}_{i}$ is differentiable and convex with respect to $\mathbf{w}$ anywhere.

(B2)

$\nabla\widetilde{l}_{i}(\mathbf{w};\mathbf{w})=\nabla l_{i}(\mathbf{w})$ for any choice of $\mathbf{w}$.

(B3)

$\widetilde{l}_{i}$ has Lipschitz continuous gradient for some constant $M$.

Taken together, these assumptions ensure that the convex surrogate loss $\widetilde{l}_{i}$ keeps the first-order properties of the original non-convex loss $l_{i}$, while allowing for a much simpler optimization. These conditions are also relatively general, allowing for a lot of flexibility in the design of the surrogate terms. Given the set of indexes $\mathcal{I}_{n}$ corresponding to the examples in the current mini-batch $\mathcal{B}_{n}$, the update at time $n$ is made by solving the following (strongly convex) surrogate optimization problem:

[TABLE]

where:

• (i)

$\mathbf{d}_{n}\in\R^{Q}$ is an auxiliary variable updated as a smoothed average of the gradients considered up to time $n$ (see below for the update equation);

• (ii)

$\rho_{n}$ is a time-dependent scalar weighting the information of the current mini-batch with respect to the historical information kept in $\mathbf{d}_{n}$;

• (iii)

the last term, with $\tau>0$, is a proximal component added to ensure that the optimization problem in (4) is strongly convex ($\tau$ can be set equal to zero if the surrogate $\widetilde{l}_{i}$, or the regularizer $r$, are already strongly convex).

For the moment we consider optimizing the entire vector $\mathbf{w}$ simultaneously, which is a special case of the algorithm in [12]; we relax this assumption later on in Section V by allowing for parallel updates of sub-blocks of the vector. Given $\widehat{\mathbf{w}}_{n+1}$, we update our current estimate with the following convex combination:

[TABLE]

where $\alpha_{n}$ is the iteration-dependent step-size. Finally, we update the auxiliary variable $\mathbf{d}_{n}$ using a similar step:

[TABLE]

where $\rho_{n}$ is the same scalar value used in the definition of (4). This update ensures that the variable asymptotically converges to the gradient of the loss term in (2). The overall algorithm is summarized in Algorithm I. Convergence to a stationary point of (2) is analyzed in the following proposition.

Proposition 1**.**

Given assumptions A-B, assume that the step-size/mixing sequences are chosen such that:

[TABLE]

Additionally, assume that the sequence $\left\{\mathbf{w}_{n}\right\}_{n=1}^{\infty}$ is bounded.222Note that this condition is trivially satisfied by imposing a finite (but arbitrarily large) box constraint guaranteeing the boundedness of the sequence. Then, all the conditions in [12, Theorem 1] are satisfied, and for every limit point generated by Algorithm I, there exists a subsequence converging to a stationary point of (2) almost surely.

Condition (iii) above ensures convergence of the auxiliary variable to the real gradient. Note that this type of almost sure convergence in stochastic non-convex settings is relatively rare in the optimization literature, which is an additional benefit of using SCA techniques [12]. To apply Algorithm I, we need a principled way to construct $\widetilde{l}_{i}(\cdot;\cdot)$. This is the topic of the next section.

IV Proposed optimization approach

IV-A Definition of the surrogate function

An immediate way to satisfy Assumption B is to define $\widetilde{l}_{i}$ as the first-order linearization of $l_{i}$:

[TABLE]

By discarding everything except first-order information, the resulting formulation does not have any definite advantage with respect to a variant of SGD. As a concrete example, consider the case of $\ell_{2}$ regularization $r(\mathbf{w})=\frac{1}{2}\left\lVert\mathbf{w}\right\rVert_{2}^{2}$. Since $r(\cdot)$ is strongly convex, we can set $\tau=0$ in (4), and we obtain the closed-form solution:

[TABLE]

where $\lambda$ is the user-defined regularization factor from (2). The resulting update resembles a simplified version of Adam [3] that does not take the second-order moment into account. A similar formulation arises when considering $\ell_{1}$ regularization instead of $\ell_{2}$. In this case, the optimum of the surrogate problem can be expressed in closed-form with the use of a soft thresholding operator, e.g., see [24, Eq. (26)].

The surrogate function in (7) destroys any information about convexity hidden in the cost function. We can do something smarter by noting that each term $l(y_{i},f(\mathbf{w};\mathbf{x}_{i}))$ in (2) is given by the composition of a non-convex function, i.e., $f(\mathbf{w};\mathbf{x}_{i})$, the NN model, with a convex loss function, i.e., $l(\cdot,\cdot)$. To preserve the convexity of the latter, we take the first-order-linearization of the NN on a single point as:

[TABLE]

where $\mathbf{J}_{i,n}=\nabla f(\mathbf{w}_{n};\mathbf{x}_{i})$ is a $Q$-dimensional vector containing the derivatives of the $i$th NN output with respect to the current estimate of the weights. More in general, it will be a matrix with one column per NN output. We refer to this quantity as the weight Jacobian. It is possible to compute it efficiently for the entire mini-batch via a single back-propagation step whose complexity is linear in the number of parameters, e.g., see [35, Section 5.3.4]. Our surrogate loss is then defined by combining $\widetilde{f}_{i}$ with the loss function, as:

[TABLE]

It is straightforward to show that (10) satisfies Assumptions B1-B3, i.e., $\widetilde{l}_{i}(\mathbf{w};\mathbf{w}_{n})$ is a differentiable, convex function that preserves the first order properties of $l_{i}$ at point $\mathbf{w}_{n}$. In particular, convexity follows from the fact that (10) is given by the composition of a convex function, i.e., $l$, with an affine mapping, i.e., $\widetilde{f}_{i}$, see [9, Section 3.2.2].

In the remainder of the section, we consider some practical examples resulting from specific choices of $l$ and $r$.

IV-B Example 1: squared loss with $\ell_{2}$ regularization

The first practical implementation we discuss is the use of a squared loss function, coupled with an $\ell_{2}$ regularization:

[TABLE]

We call this the ridge regression cost in analogy with the linear case. This is the most common way of training neural networks for regression, e.g., see [36] for a very recent example. As before, we can set $\tau=0$ thanks to the presence of a strongly convex regularization term. We define for convenience the ‘residual’ terms $r_{i,n}$ as:

[TABLE]

After simple algebra manipulations, we can write the optimization problem in (4) as a quadratic optimization problem:

[TABLE]

with $\mathbf{I}$ being the identity matrix of appropriate size, and we defined:

[TABLE]

Thus, the solution of the quadratic problem (13) is given by:

[TABLE]

This requires the inversion of a $Q\times Q$ matrix, which can become impractical for large $Q$. In Section V we show a principled way to decompose this problem and obtain a significant speedup when using multiple processors. Alternatively, one can solve the original problem in (13) using, e.g., highly customized conjugate gradient optimization procedures [37]. If we set $\rho_{n}=1$ and $\mathcal{B}_{n}=\mathcal{S}$ for any $n$ (i.e., we work in a batch fashion), we recover the PL-SCA algorithm proposed in [24], where PL stands for ‘partial linearization’. For a very small $\lambda$, solving (16) can give rise to numerical problems, because $\mathbf{A}_{n}$ has at most rank $L$, being a sum of $L$ rank-1 matrices. In this case, one can set $\tau>0$, obtaining a slightly modified solution:

[TABLE]

which still guarantees convergence of the stochastic SCA procedure to a local solution of (2).

IV-C Case 2: sparsity-inducing penalties

As a second use case, we consider again the use of a squared loss function, this time combined with an $\ell_{1}$ regularization to achieve sparsity on the weights:

[TABLE]

Proceeding as in (13), we can immediately formulate the surrogate problem as an $\ell_{1}$-regularized quadratic problem:

[TABLE]

There is a wide range of fast solvers for (19) [38], most notably the fast iterative shrinkage and thresholding algorithm (FISTA) [39], which achieves a $\mathcal{O}(1/n^{2})$ rate of convergence. Note that this approach will yield exactly sparse solutions at every iteration. On the contrary, it is customary in the NN literature to solve $\ell_{1}$ regularized problems with SGD algorithms, in which case an exactly sparse solution can never be reached [40], and a further thresholding step is needed.

This formulation can be extended immediately to similar forms of $\ell_{1}$ regularization, such as elastic net penalties (given by a weighted sum of $\ell_{2}$ and $\ell_{1}$ normalization), and recent group sparse formulations such as [41] and [42]. Group sparse penalties can be used to favor structured forms of sparsity where entire neurons are removed in the optimization process. Suppose the neurons are indexed as $1,\ldots,P$, and denote by $\mathbf{w}_{p}\subset\mathbf{w}$ the set of weights outgoing from the $p$th neuron, such that:

[TABLE]

Group sparse regularization is achieved as:

[TABLE]

where $a_{p}$ are scalar coefficients defined as the square root of the dimensionality of the corresponding groups. Most optimization algorithms designed to optimize (19) can be applied equivalently even if we interchange the $\ell_{1}$ term with the group sparse term in (21). Again, this allows us to obtain exactly sparse solutions.

IV-D Case 3: cross-entropy loss

As a third use case, we consider a binary classification problem with $y_{i}=\left\{0,1\right\}$, for which it is common to optimize the cross-entropy loss defined as:

[TABLE]

The previous loss can be combined with either $\ell_{2}$ or $\ell_{1}$ regularization, depending on the learning task. Some care must be taken here because, even when the original NN model $f(\cdot)$ is always bounded, the same is not true for its linearization in (9). Substituting (9) in (22) is then undefined whenever $\widetilde{f}$ is non-positive or larger than $1$. To solve this issue, note that for binary classification problems with the cross-entropy loss the NN can always be written as:

[TABLE]

where $\sigma(\cdot)$ is a squashing function (e.g., sigmoid) ensuring that the output is properly defined as a probability, and $f^{L}$ denotes the output of the NN up to the last nonlinearity. We define a linearization on $f^{L}$ similar to (9):

[TABLE]

where $\mathbf{J}^{L}_{i,n}=\nabla f^{L}(\mathbf{w}_{n};\mathbf{x}_{i})$. The sigmoid function is neither convex nor concave, but the combination of a sigmoid with (22) is convex. Thus, proceeding as in (10), a proper surrogate to use in (4) is obtained as:

[TABLE]

The final optimization problem in (4) is similar to a standard logistic regression, with the straightforward inclusion of a linear term and a proximal norm on $\mathbf{w}$. Also here we can exploit highly customized solvers for optimizing the resulting strongly convex problem, e.g., see [43] for a recent survey up to $2012$.

IV-E Case 4: non-convex regularizers

The previous sections described use cases that are relatively common in the NN literature. One interesting extension concerns the use of non-convex regularization terms. These terms generally arise because of the need to regularize (in some meaningful way) the output of the NN, or the activations of specific neurons. By making both the loss term and the regularization term highly non-convex, they are generally harder to optimize and less common in the literature. However, they fit naturally in the SCA framework, because we can apply the same ideas described in Section IV to also convexify the regularization term.

As a specific case, consider the use of manifold regularization in deep networks [44]. The basic idea is to force the NN to provide similar outputs whenever two inputs are ‘close’ according to some distance measure. To this end, suppose that $q_{ij}$ is a non-negative value measuring the distance between the inputs $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$. Typically, this is defined as some measure of the Euclidean distance for the $k$-nearest neighbors of $\mathbf{x}_{i}$, and [math] otherwise [44]. A manifold regularization term can then be written as:

[TABLE]

Note that, for a mini-batch of elements, computation of (26) requires $Lk$ additional forward/backward computations in general [44], one for each neighbor of the elements in the mini-batch. We can handle this sort of regularization by replacing $r(\cdot)$ in (4) with a strongly convex approximation, following the same methodology as in Section IV-A. In particular, substituting $f(\cdot)$ with its first-order linearization in (26) we can rewrite it as:

[TABLE]

where we defined:

[TABLE]

Thus, manifold regularization provides just a linear and a quadratic term in $\mathbf{w}$, and can easily be plugged-in with any choice of loss function described before. Particularly, a closed-form solution is preserved when employing the squared loss, whereas a simple addition of a quadratic term is achieved in the case of the cross-entropy loss function.

V Parallelizing the surrogate optimization

In general, the algorithms described in the previous section are more expensive in computational terms with respect to accelerated gradient methods. This is both their advantage (resulting in faster convergence speed), and their drawback when moving to a large-scale regime. In this section, we show a simple way to parallelize their computation, whenever we have access to a multi-core (or multi-machine) environment.

Roughly, NN training algorithms can be parallelized by partitioning the data (i.e., providing a different mini-batch to each processing unit), or by partitioning the optimization variable $\mathbf{w}$. For the moment, we focus on the latter strategy. To this end, suppose that $\mathbf{w}$ is partitioned in $C$ non-overlapping blocks $\mathbf{w}_{1},\ldots,\mathbf{w}_{C}$, so that $\mathbf{w}=\bigcup_{c=1}^{C}\mathbf{w}_{c}$. $\mathbf{w}_{-c}\triangleq(\mathbf{w}_{p})_{1=p\neq c}^{C}$ will denote the tuple of all blocks excepts the $c$-th one, and similarly for all other variables. Additionally, we assume that the regularization term $r$ is block separable, i.e., $r(\mathbf{w})=\sum_{c=1}^{C}r_{c}(\mathbf{w}_{c})$ for some $r_{c}$. This is true for the $\ell_{2}$ and $\ell_{1}$ norms, and it holds true also for the group sparse norm in (21) if we choose the groups in a consistent way.

At the $n$th iteration each computing unit gets assigned the $c$th block $\mathbf{w}_{c,n}$ of $\mathbf{w}_{n}$. Then, each core solves a smaller surrogate problem defined as:

[TABLE]

where $\widetilde{l}_{i,c}(\cdot)$ is a surrogate term respecting assumption B on the block $\mathbf{w}_{c}$ only. Each core $c$ can then minimize its corresponding term independently of the others, and their solutions can be aggregated to form the final solution vector. To obtain the surrogate function associated to each core $c$, we simply compute a full surrogate as in the previous section, then fix all the variables $\mathbf{w}_{-c,n}$ to their current value, such that the resulting function depends only on $\mathbf{w}_{c}$.

As an example of this idea, consider the ridge surrogate defined in Section IV-B. Simple algebra shows that the solution of the local surrogates in a parallel environments are given by:

[TABLE]

where $\mathbf{A}_{c,c,n}$ is the block (rows and columns) of the matrix $\mathbf{A}_{n}$ in (14) corresponding to the $c$-th partition, whereas $\mathbf{A}_{c,-c,n}$ takes the rows corresponding to the $c$-th partition and all the columns not associated to $c$. Each core has now to invert a matrix having (approximately) size $\frac{1}{C}$ of the original one, thus remarkably reducing the overall computational burden. Similar arguments can be used also to parallelize all the other formulations described in the previous section.

The theorems in [12] ensure that convergence is guaranteed even in the parallel case. Interestingly, convergence is also guaranteed in the more general case where the number of blocks is larger than the number of computational units, and at each iteration every processor is randomly assigned a block [12]. Thus, the size of the blocks (and, consequently, of the surrogate problems to be solved) can be freely chosen based on the available computational requirements. Even more generally, one can consider randomly assigning both a block variable and a separate mini-batch. Note, however, that in general every core will have to compute the entire Jacobian matrix in (9) due to the back-propagation step. As a matter of fact, this is a general limitation of any parallel gradient-based approach to NN training [45], which is generating new fields of research in term of signal back-propagation, e.g., see [46].

VI Experimental results

In this section we evaluate the convergence behavior of the proposed technique. We focus on the ridge case considered in Section IV-B. We consider several medium-sized datasets in Section VI-A, and a large-scale problem in Section VI-B. Python code to repeat the experiments is available under open-source license on the web.333https://bitbucket.org/ispamm/sca-optimization-for-neural-networks Back-propagation and part of the comparing algorithms are based on the AutoGrad library [47].

For all experiments, the original dataset is normalized so that inputs lie in the $\left[-0.5,0.5\right]$ range, and outputs lie in the $\left[-0.9,0.9\right]$ range. For each run, a random $25\%$ of the dataset is used for testing, and the rest for training the network. All experiments are repeated $100$ times by varying the data partitioning and the NN initialization. Missing data is replaced with the median value for the corresponding feature in the entire dataset.

Regarding the NN structure, we use hyperbolic tangent nonlinearities in all neurons, and weights are initialized using the normalized strategy described by [48].

VI-A Experiments on mid-sized datasets

We start by consider four mid-sized regression datasets, whose characteristics are briefly summarized in Table II. All of them were downloaded from the UCI repository.444https://archive.ics.uci.edu/ml/ The fourth column in Table II describe the topology of the NN we have chosen. These parameters are chosen based on an analysis of previous literature in order to obtain state-of-the-art results. However, we underline that our aim is to compare different solvers for the same NN optimization problem, and for this reason only relative differences in accuracy are of concern. For all datasets, we select a small regularization coefficient $\lambda=10^{-3}$, which is found to provide good results.

We compare the results of the algorithm in Section IV-B with respect to four state-of-the-art solvers, in terms of mean-squared error (MSE) over the test data, when solving the global optimization problem with the ridge regression cost in (11). Specifically, we consider the following algorithms:

Stochastic gradient descent (GD)

: this is a simple first-order steepest descent procedure with diminishing step-size (see below).

Adagrad

: differently from SGD, we use different step-sizes per weight, which evolve according to the relative values of the gradients’ updates [49].

RMSProp

: it also considers adaptive independent step-sizes; however they are adapted based on an exponentially-weighted moving average [50].

Adam

: Adam combines a momentum strategy with adaptive step-sizes, where both first- and second-order moments are computed in an streaming fashion [3].

For the proposed algorithm and SGD, we use the following quadratically decreasing rule for selecting the step-sizes:

[TABLE]

where $\alpha_{0}$ and $\varepsilon$ are selected by the user. The same sequence is also used for selecting $\rho_{n}$ in the proposed algorithm. For fairness of comparison, we selected default values for all algorithms leading to (in average) their best convergence behavior. In particular, we set $\alpha_{0}=0.5$, $\varepsilon=0.01$, $\tau=0$, and $\rho_{0}=0.9$ for SCA, $\alpha_{0}=0.1$, and $\varepsilon=0.01$ for SGD. For AdaGrad and RMSProp, we set the initial learning rate to $0.01$. For the latter, we set a decaying value of $\gamma=0.9$. For Adam, we use the default values as in [3]. For all algorithms, we consider randomly extracted mini-batches of size $L=20$.

Results in terms of overall loss value per iteration are shown in Fig. 1, where the proposed algorithm is shown with a purple line, and the shaded areas represent the standard deviation around the mean. It can be seen that the overall results are generally consistent across the four datasets. Specifically, SGD is a relatively poor choice, getting stuck in a close minimum two out of four times (CASP and Wine datasets), and converging slowly in the other two cases (Parkinson and Skills). Among the state-of-the-art algorithms, Adam is always the best one (which is in accordance with its popularity), closely followed by Adagrad, with RMSProp scoring in the middle between these two and SGD. However, in all situations the proposed SCA technique is able to significantly outperform the other algorithms, being consistently faster in two out of four cases (Skills and Wine), and converging to a better solution in all settings.

The better results in term of convergence are equivalently found when considering the MSE on the independent test set, as shown in Table III. It can be seen that the final MSE for the SCA training algorithm always outperforms the MSE for the architectures optimized by competing algorithms. Considering the training time, we note that for these algorithms all algorithms required (approximately) the same training time per iteration, which is not shown here for reasons of space.

An interesting question is motivating theoretically the improvement in convergence time provided by the SCA methodology. To this end, one can observe that the matrix (14) is in fact an approximation to the true Hessian matrix of the squared cost, which is obtained by assuming that the error is uncorrelated with the second derivative of the squared loss function. This is known as an outer-product approximation, or Levenberg-Marquardt approximation (see [35, Section 5.4.2] for a discussion). Thus, it is reasonable to assume that the overall algorithm is able to maintain some information on the curvature of the cost function, even if higher-order derivatives with respect to the gradient are never explicitly computed. A similar argument was made in [24]. This is an interesting line of reasoning which could eventually lead to improved approximations for the cost function.

VI-B Experiment on a large-scale dataset

Next, we consider a large-scale dataset to evaluate the performance of the algorithm on high-dimensional NNs. Due to the size of the parameter vector, this also allows us to test the parallel version of the algorithm discussed in Section V. To this end, we consider the SUSY benchmark originally introduced in [51].555https://archive.ics.uci.edu/ml/datasets/SUSY The dataset is composed by $5$ million simulations of particle collisions at high energy, simulating the environment found in currently used particle accelerators. The task is to distinguish between processes where supersymmetric particles are created (denoted as $\chi^{\pm}$ and $\chi^{0}$), and background processes. The challenge is that supersymmetric particles are not observed, and in both cases the observed particles are identical (leptons). Each example is described by an $18$-dimensional input vector, where the first $8$ features are low-level features describing the measurements, while the remaining $10$ features are high-level features constructed from the low-level ones. For this experiment, we use larger mini-batches of $L=50$ elements. The network has two hidden layers with $100$ neurons each, and we use a slightly larger regularization factor $\lambda=0.01$.

In Fig. 2 we show the convergence per iteration of Adam as compared to the proposed approach. The initial learning rates for the two algorithms were fine-tuned to obtain the fastest convergence behavior. Also, we use a parallel version with $C=4$, whose training time is roughly twice with respect to Adam, which is acceptable. It can be seen that, as for the previous section, the SCA technique significantly outperforms the Adam approach. The average test accuracy is similarly larger: on an independently kept $25\%$ of the original dataset, the network trained via the SCA algorithm has an average area under the curve (AUC) of $0.85$, which is eight points higher than the one obtained by Adam ($0.77$).

To visualize the speedup obtained by the parallelization procedure, in Fig. 3 we show the relative speedup with respect to $C=1$ obtained when varying $C$ in $2,4,8,\ldots,64$. As can be shown, the training time for $C=2$ is generally higher than Adam, but the gap closes substantially for $C=4$. The two algorithms have almost the same training time for $C=8$, while SCA requires less than half the training time of Adam when having available $C=64$ computational units. This speedup is obtained without sacrificing accuracy. In Fig. 4, we show the receiver operating characteristic (ROC) curves obtained by Adam, SCA with $C=4$, and SCA with $C=16$. We can see that SCA outperforms the Adam algorithm (in line with our previous discussion). Additionally, the same ROC curve is obtained for both configurations of the method.

VII Conclusions

In this paper, we proposed a novel family of algorithms designed for optimizing neural networks. At each iteration of the algorithm, a random mini-batch of data is extracted from the training set, and a strongly convex approximation of the original training algorithm is solved. The algorithm only requires first-order information on the network. Our experimental results show that it performs favorably with respect to state-of-the-art approaches, being in general faster to converge to a better minimum of the optimization problem. For large-scale problems, the algorithm can be easily parallelized across multiple computational units. Further research will investigate the possibility of designing better approximation functions, and the customization of the framework to different families of NNs, including convolutional and recurrent networks. Additionally, we plan to test the algorithm on a large-scale cluster environment with asynchronous updates.