Stochastic Trust Region Inexact Newton Method for Large-scale Machine Learning

TL;DR

This paper introduces STRON, a stochastic trust region inexact Newton method that combines subsampling, variance reduction, and conjugate gradient techniques to efficiently solve large-scale machine learning problems.

Contribution

The paper proposes a novel stochastic second-order optimization method, STRON, integrating subsampling and trust region strategies for large-scale learning.

Findings

STRON outperforms existing methods on benchmark datasets.

Variance reduction techniques did not improve large-scale learning.

Preconditioned conjugate gradient was less effective in this context.

Abstract

Nowadays stochastic approximation methods are one of the major research direction to deal with the large-scale machine learning problems. From stochastic first order methods, now the focus is shifting to stochastic second order methods due to their faster convergence and availability of computing resources. In this paper, we have proposed a novel Stochastic Trust RegiOn Inexact Newton method, called as STRON, to solve large-scale learning problems which uses conjugate gradient (CG) to inexactly solve trust region subproblem. The method uses progressive subsampling in the calculation of gradient and Hessian values to take the advantage of both, stochastic and full-batch regimes. We have extended STRON using existing variance reduction techniques to deal with the noisy gradients and using preconditioned conjugate gradient (PCG) as subproblem solver, and empirically proved that they do not work as expected, for the large-scale learning problems. Finally, our empirical results prove efficacy of the proposed method against existing methods with bench marked datasets.

Tables

Table 1. Table 1: Datasets used in experimentation

Dataset	#features	#datapoints
gisette	5000	6,000
rcv1.binary	47,236	20,242
webspam (unigram)	254	350,000
covtype.binary	54	581,012
ijcnn1	22	49,990
news20.binary	1,355,191	19,996
real-sim	20,958	72,309
Adult	123	32,561
mushroom	112	8124

Table 2. Table 2: Comparison of Training Time (seconds) for low accuracy ( ϵ italic-ϵ \epsilon =0.01) solution

Datasets$\downarrow$	Methods$\to$	SVRG_LBFGS	STRON	SVRG_SQN	Newton_CG	TRON
covtype	Time	0.8440$\pm$0.1770	0.0620$\pm$0.0045	0.4140$\pm$0.0397	0.9740$\pm$0.1907	0.7294$\pm$0.0162
	Accuracy	1.00$\pm$0.0	1.00$\pm$0.0	1.00$\pm$0.0	1.00$\pm$0.0	1.00$\pm$0.0
real-sim	Time	1.5920$\pm$0.5867	0.9340$\pm$0.1258	2.3480$\pm$0.5909	0.68$\pm$0.01	0.88$\pm$0.10
	Accuracy	0.9670$\pm$0.0055	0.9697$\pm$0.0017	0.9670 $\pm$0.0029	0.9688$\pm$0.0022	0.9698$\pm$0.0015
rcv1	Time	2.1780$\pm$0.1813	0.5700$\pm$0.0200	3.2840$\pm$0.0270	0.5860$\pm$0.0167	0.6589$\pm$0.0866
	Accuracy	0.9641$\pm$0.0023	0.9640$\pm$0.0025	0.9637$\pm$0.0024	0.9640$\pm$0.0024	0.9639$\pm$0.0027
news20	Time	42.8960$\pm$2.5993	6.9620$\pm$1.3534	84.0520$\pm$2.2915	6.764$\pm$0.0074	7.2470$\pm$0.3567
	Accuracy	0.9333$\pm$0.0079	0.9337$\pm$0.0076	0.9339$\pm$0.0072	0.9338$\pm$0.0074	0.9338$\pm$0.0077
mushroom	Time	0.1540$\pm$0.01340	0.0820$\pm$0.0268	0.1820 $\pm$0.01920	0.072$\pm$0.0045	0.1316$\pm$0.0169
	Accuracy	0.9995$\pm$0.0007	0.9992$\pm$0.0005	0.9986$\pm$0.0026	0.9996$\pm$0.0005	0.9995$\pm$0.0005
ijcnn1	Time	0.3460$\pm$0.1390	0.2560$\pm$0.03780	0.3600$\pm$0.1158	0.3120$\pm$0.0795	0.2564$\pm$0.0185
	Accuracy	0.9235$\pm$0.0022	0.9232$\pm$0.0021	0.9237$\pm$0.0018	0.9238$\pm$0.0025	0.9233$\pm$0.0185
Adult	Time	0.3620$\pm$0.0981	0.2000$\pm$0.0187	0.3100$\pm$0.0644	0.2380$\pm$0.0084	0.2139$\pm$0.0046
	Accuracy	0.8479$\pm$0.0028	0.8470$\pm$0.0028	0.8469$\pm$0.0031	0.8473$\pm$0.0030	0.8477$\pm$0.0028
gisette	Time	9.5620$\pm$0.9891	8.1800$\pm$0.4260	10.7020$\pm$1.1314	41.6480$\pm$36.1190	10.4277$\pm$0.1709
	Accuracy	0.9725$\pm$0.0036	0.9730$\pm$0.0045	0.9690$\pm$0.0026	0.9701$\pm$0.0057	0.9755$\pm$0.0028
webspam	Time	6.0480$\pm$2.2814	6.3720$\pm$0.9942	7.6080$\pm$2.1462	3.272$\pm$0.2141	9.0548$\pm$1.2572
	Accuracy	0.9139$\pm$0.0033	0.9214$\pm$0.0019	0.9158$\pm$0.0036	0.9219$\pm$0.0017	0.9227$\pm$0.0009

Equations

$$\underset{w}{\min }R(w)=\mathbb{E}\left[f\left(w;\xi \right)\right],$$

$$\underset{w}{\min }F(w)=\frac{1}{l}\sum _{i=1}^{l}f(w;x_i,y_i).$$

$$w^{k+1}=w^k-{\alpha }_k\nabla F\left(w_k\right),$$

$$w^{k+1}=w^k-{\alpha }_k\nabla f_{i_k}\left(w_k\right),$$

$$w^{k+1}=w^k-{\alpha }_k{\nabla }^{2}F{\left(w_k\right)}^{-1}\nabla F\left(w_k\right).$$

$$w^{k+1}=w^k-{\alpha }_k{B}_k\nabla F\left(w_k\right),$$

$$\begin{array}{r}F(w_k+p)-F(w_k)\approx m_k(p)\equiv \nabla F{\left(w_k\right)}^{T}p+\frac{1}{2}{p}^T{\nabla }^{2}F\left(w_k\right)p.\end{array}$$

$$\begin{array}{r}{\nabla }^{2}F\left(w_k\right)p=-\nabla F\left(w_k\right),\end{array}$$

$$\begin{array}{r}{w}^{k+1}=w^k+p_k=w^k-{\nabla }^{2}F{\left(w_k\right)}^{-1}\nabla F\left(w_k\right).\end{array}$$

$$\Vert r_k\Vert \le {\eta }_k^{{}^{\prime }}\Vert \nabla F\left(w_k\right)\Vert ,$$

$$\underset{p}{\min }{m}_k\left(p\right)\text{s.t.}\Vert p\Vert \le {△}_k,$$

$${\rho }_k=\frac{F(w_k+p_k)-F(w_k)}{m_k(p_k)}.$$

$$w_{k+1}=\{\begin{array}{ll}{w}_k+p_k,& \text{if}{\rho }_k>{\eta }_0,\\ w_k,& \text{if}{\rho }_k\le {\eta }_0,\end{array}$$

$${△}_{k+1}\in \{\begin{array}{ll}\left[{\gamma }_1\min \{\Vert p_k\Vert ,{△}_k\},{\gamma }_2{△}_k\right],& \text{if}{\rho }_k\le {\eta }_1,\\ \left[{\gamma }_1{△}_k,{\gamma }_3{△}_k\right],& \text{if}{\rho }_k\in \left({\eta }_1,{\eta }_2\right),\\ \left[{△}_k,{\gamma }_3{△}_k\right],& \text{if}{\rho }_k\ge {\eta }_2,\end{array}$$

$$\begin{array}{r}\underset{p}{\min }{m}_k(p)=\nabla F_{X_k}{\left(w_k\right)}^{T}p+\frac{1}{2}{p}^T{\nabla }^2{F}_{S_k}\left(w_k\right)p,\text{s.t.}\Vert p\Vert \le {△}_k,\end{array}$$

$$\begin{array}{r}{\nabla }^2{F}_{S_k}\left(w_k\right)=\frac{1}{|S_k|}\sum _{i\in S_k}{\nabla }^2{f}_i\left(w_k\right),\\ \nabla F_{X_k}\left(w_k\right)=\frac{1}{|X_k|}\sum _{i\in X_k}\nabla f_i\left(w_k\right),\\ F_{X_k}\left(w_k\right)=\frac{1}{|X_k|}\sum _{i\in X_k}{f}_i\left(w_k\right),\end{array}$$

$$\Vert r_k\Vert \le {\eta }_k^{{}^{\prime }}\Vert \nabla F_{X_k}\left(w_k\right)\Vert ,$$

$$\underset{w}{\min }F(w)=\frac{1}{l}\sum _{i=1}^l\log \left(1+\exp \left(-y_i{w}^T{x}_i\right)\right)+\frac{\lambda }{2}\Vert w{\Vert }^2.$$

$$\underset{w}{\min }F(w)=\frac{1}{l}\sum _{i=1}^l\max (0,1-y_i{w}^T{x}_i{)}^2+\frac{\lambda }{2}\Vert w{\Vert }^2.$$

$$M^{-1}{\nabla }^{2}F(w)p=-M^{-1}\nabla F(w).$$

$$M=\alpha \times diag(H)+(1-\alpha )\times I,$$

$$\underset{\hat{p}}{\min }{\left(L^{-1}\nabla F_{X_k}\left(w_k\right)\right)}^T\hat{p}+\frac{1}{2}{\hat{p}}^T\left(L^{-1}{\nabla }^2{F}_{S_k}\left(w_k\right)L^{-T}\right)\hat{p},\text{s.t.}\Vert \hat{p}\Vert \le {△}_k,$$

$$\begin{array}[]{ll}&O\left(n\right)\times\#\text{CG iterations}+O(nl).\end{array}$$

$$\begin{array}[]{ll}&O\left(n\right)\times\#\text{CG iterations}+O(n|S_{k}|).\end{array}$$

$$g_k=\nabla F_{S_k}(w_k)-\nabla F_{S_k}(\bar{w})+\nabla F(\bar{w}),$$

Full text

1 Introduction

Machine learning involves data intensive optimization problems which have large number of component functions corresponding to large amount of available data. Traditional/classical methods, like Newton method, fail to perform well on such large-scale learning problems (i.e., problems with large number of data points) due to large per-iteration complexity. So nowadays one of the major challenge in machine learning is to develop scalable and efficient algorithms to deal with these large-scale learning problems Zhou et al. (2017); Chauhan et al. (2017, 2018a).

To solve the machine learning problems, gradient descent (GD) Cauchy (1847) is the classical method of choice but it trains slowly while dealing with large-scale learning problems due to high per-iteration cost. Stochastic approximation based methods Robbins and Monro (1951) can be quite effective in such situations but they converge slowly due to the noisy approximations of gradient. So a variety of stochastic variance reduction techniques came to existence, e.g., Le Roux et al. (2012); Defazio et al. (2014); Schmidt et al. (2016); Allen-Zhu (2017); Fanhua et al. (2018). But the major limitation of these methods is that they can converge up to linear rate only.

Newton method is another classical method to solve optimization problems, which can give up to quadratic convergence rate Boyd and Vandenberghe (2004). But again (pure) Newton method is not feasible with large-scale learning problems due to huge per-iteration computational complexity and need to store a huge Hessian matrix. So nowadays one of the most significant open question in optimization for machine learning is: “Can we develop stochastic second order methods with quadratic convergence, like Newton method but has low per-iteration complexity, like stochastic approximation methods?”. After success in the stochastic first order methods, the research is shifting its focus towards the stochastic second order methods to leverage the faster convergence of second order methods and the available computing power.

Inexact Newton (also called truncated Newton or Hessian free) methods and quasi-Newton methods are among the major research directions for developing second order methods Bollapragada et al. (2016). Inexact Newton methods try to solve the Newton equation approximately without calculating and storing the Hessian matrix. On the other hand, quasi-Newton methods try to approximate the Hessian inverse and avoid the need to store the Hessian matrix. Thus both the methods try to resolve the issues with Newton method for solving large-scale learning problems. The stochastic variants of inexact Newton and quasi-Newton, further reduce the complexity of these methods by using subsampled gradient and Hessian calculations. In this paper, we have proposed a novel stochastic trust region inexact Newton (STRON) method to solve the large-scale learning problems, which introduces subsampling to gradient and Hessian calculations. It uses progressive subsampling to enjoy the benefits of both the regimes, stochastic approximation and full-batch learning. We further extend the method using existing variance reduction techniques to deal with noise produced by subsampling of gradient values, and by proposing PCG for solving the trust region subproblem.

1.1 Optimization Problem

We consider unconstrained convex optimization problem of expected risk, as given below:

[TABLE]

where $f\left(w;\xi\right)=f\left(w;x_{i},y_{i}\right)=f\left(h\left(w;x_{i}\right),y_{i}\right)$ is a smooth composition of linear prediction model $h$ and loss function $f$ over randomly selected data point $\left(x_{i},y_{i}\right)$ from the unknown distribution $P\left(x_{i},y_{i}\right)$, parameterized by the model parameter $w\in\mathbb{R}^{n}$. Since it is not feasible to solve (1) as $P$ is unknown so the model is approximated by taking a set $N=\left\{(x_{1},y_{1}),...,(x_{l},y_{l})\right\}$ of $l$ data points from the unknown distribution $P$ and then solving the empirical risk minimization problem, as given below:

[TABLE]

For simplicity, we take $f\left(w;x_{i},y_{i}\right)=f_{i}(w)$. Finite sum optimization problems of type (2) exists across different fields, like signal processing, statistics, operation research, data science and machine learning, e.g., logistic regression and SVM in machine learning.

1.2 Solution Techniques

Simple iterative classical method to solve (2) is gradient descent (GD) Cauchy (1847), as given below:

[TABLE]

where $(k+1)$ is iteration number and $\alpha_{k}$ is called learning rate or step size. The complexity of this iteration is $O\left(nl\right)$ which is large for large-scale learning problems due to large number of data points. Stochastic gradient descent (SGD) Robbins and Monro (1951) is very effective to deal with such problems due to its low per-iteration complexity, as given below:

[TABLE]

where $i_{k}$ is randomly selected data point. But convergence can’t be guaranteed in SGD because of noisy gradient values.

Another classical second order method to solve (2) is Newton method, as given below:

[TABLE]

The complexity of this iteration is $O(n^{2}l+n^{3})$ and it needs to store and invert the Hessian matrix $\nabla^{2}F\left(w_{k}\right)$, which is computationally very expensive and needs large memory for large-scale learning problems, respectively. That’s why first order methods and their stochastic variants have been studied very extensively, during the last decade, to solve the large-scale learning problems but not second order methods. As stochastic first order methods have hit their limits and due to good availability of computing power, the main focus is shifting towards stochastic second order methods, and nowadays one important open question is to find stochastic second order methods with quadratic convergence rates.

There are two major research directions for second order methods: quasi-Newton methods and inexact Newton methods, both of which try to resolve the issues associated with the Newton method. Quasi-Newton methods try to approximate the Hessian matrix during each iteration, as given below:

[TABLE]

where $B_{k}$ is an approximate of $\nabla^{2}F\left(w_{k}\right)^{-1}$, e.g., Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm is one such method Fletcher (1980). On the other hand, inexact Newton methods try to solve the Newton system approximately, e.g., Newton-CG Steihaug (1983). Both the methods try to resolve the issues related with Newton method but still their complexities are large for large-scale problems. So a number of stochastic variants of these methods have been proposed, e.g., Byrd et al. (2011, 2016); Bollapragada et al. (2018); Bellavia et al. (2018) which introduce subsampling to gradient and Hessian calculations.

1.3 Contributions

The contributions of the manuscript are listed below:

• •

The manuscript highlights the recent shift of stochastic methods from first order to second order methods and raises an open question that we need stochastic second order methods with quadratic convergence rate to solve the large-scale machine learning problems.

• •

We have proposed a novel subsampled variant of trust region inexact Newton method for solving large-scale learning problems, which is called STRON and is the first stochastic variant of trust region inexact Newton methods. STRON uses progressive subsampling scheme for gradient and Hessian calculations to enjoy the benefits of both stochastic and full batch regimes. STRON can converge up to quadratic rate and answers the raised open question.

• •

STRON has been extended using existing variance reduction techniques to deal with the noisy approximations of the gradient calculations. The extended method uses stochastic variance reduced gradient (SVRG) for variance reduction with static batching for gradient calculations and progressive batching for the Hessian calculations. The empirical results prove that this does not work as expected for large-scale learning.

• •

We further extend STRON and use PCG, instead of CG method, to solve the Newton system inexactly. We have used weighted average of identity matrix and diagonal matrix as the preconditioner. But even this fails to work as expected for large-scale learning.

• •

Finally, our empirical experiments prove the efficacy of STRON against existing techniques with bench marked datasets.

2 Literature Review

Stochastic approximation methods are very effective to deal with the large-scale learning problems due to their low per-iteration cost, e.g., SGD Robbins and Monro (1951), but they lead to slow convergence rates due to the noisy approximation. To deal with the noise issue, a number of techniques have been proposed and some of the important techniques (as discussed in Csiba and Richt (2016)) are: (a) decreasing learning rates Shalev-Shwartz et al. (2007), (b) using mini-batching Chauhan et al. (2018b), (c) importance sampling Csiba and Richt (2016), and (d) variance reduction Johnson and Zhang (2013). The variance reduction methods can further be classified into three categories: primal methods Schmidt et al. (2016); Chauhan et al. (2019b), dual methods Shalev-Shwartz and Zhang (2013) and primal-dual methods Zhang and Xiao (2015). The variance reduction techniques are effective to deal with the large-scale learning problems because of low per-iteration complexity, like SGD, and have fast linear convergence, like GD. These techniques exploit the best of SGD and GD but for these stochastic variants of first order methods the convergence is limited to linear rate only, unlike the second order methods which can give up to quadratic rate.

Second order methods utilize the curvature information to guide the step direction towards the solution and exhibit faster convergence than the first order methods. But huge per-iteration cost due to the huge Hessian matrix and its inversion make the training of models slow for large-scale problems. So certain techniques have been developed to deal with the issues related to Hessian matrix, e.g., quasi-Newton methods and inexact Newton methods are two major directions to deal with the huge computational cost of Newton method. Quasi-Newton methods approximate the Hessian matrix during each iteration, e.g., BFGS Fletcher (1980) and its limited memory variant, called L-BFGS Liu and Nocedal (1989), are examples of the quasi-Newton class which use gradient and parameter values from the previous iterations to approximate the Hessian inverse. L-BFGS uses only recent information from previous $M$-iterations. On the other hand, inexact Newton methods try to solve the Newton system approximately, e.g., Newton-CG Steihaug (1983).

Recently, several stochastic variants of BFGS and L-BFGS have been proposed to deal with large-scale problems. Schraudolph et al. (2007) proposed stochastic variants of BFGS and L-BFGS for the online setting, known as oBFGS. Mokhtari and Ribeiro (2014) extended oBFGS by adding regularization which enforces upper bound on the eigen values of the approximate Hessian, known as RES (Regularized Stochastic BFGS). Stochastic quasi-Newton (SQN) Byrd et al. (2016) is another stochastic variant of L-BFGS which collects curvature information at regular intervals, instead of at each iteration. Variance-reduced Stochastic Newton (VITE) Lucchi et al. (2015) extended RES and proposed to use variance reduction for the subsampled gradient values for solving smoothly strongly convex problems. Kolte et al. (2015) provided another stochastic L-BFGS method with variance reduction using SVRG (referred as SVRG-LBFGS). Moritz et al. (2016) proposed Stochastic L-BFGS (SLBFGS) using SVRG for variance reduction and using Hessian-vector product to approximate the gradient differences for calculating the Hessian approximations, also referred as SVRG-SQN. SVRG-LBFGS and SVRG-SQN differ in how the curvature is approximated (i.e. how Hessian is approximated), former collects the curvature information once during each epoch (outer iterations) using gradient differences but later collects the curvature information after regular intervals inside the inner-iterations using Hessian-vector products. Berahas et al. (2016) proposed multi-batch scheme into stochastic L-BFGS where batch sample changes with some overlaps with previous iteration. Bollapragada et al. (2018) proposed progressive batching, stochastic line search and stable Newton updates for L-BFGS. Bollapragada et al. (2016) studies the conditions on the subsample sizes to get the different convergence rates.

Stochastic inexact Newton methods are also explored extensively. Byrd et al. (2011) proposed stochastic variants of Newton-CG along with L-BFGS method. Bollapragada et al. (2016) studies subsampled Newton methods and find conditions on subsample sizes and forcing term (constant used with the residual condition), for linear convergence of Newton-CG method. Bellavia et al. (2018) studies the effect of forcing term and line search to find linear and super-linear convergence of Newton-CG method. Newton-SGI (stochastic gradient iteration) is another way of solving the linear system approximately and is studied in Agarwal et al. (2017).

Trust Region Newton (TRON) method is one of the most efficient solver for solving large-scale linear classification problems Lin et al. (2008). This is trust region inexact Newton method which does not use any subsampling and is present in LIBLINEAR library Fan et al. (2008). Hsia et al. (2017) extends TRON by improving the trust region radius value. Hsia et al. (2018), further extends TRON and uses preconditioned conjugate gradient (PCG) which uses weighted average of identity matrix and diagonal matrix as a preconditioner, to solve the trust region subproblem. Since subsampling is an effective way to deal with the large-scale problems so in this paper, we have proposed a stochastic variant of trust region inexact Newton method, which have not been studied so far to the best of our knowledge.

3 Trust Region Inexact Newton Method

Inexact Newton methods, also called as Truncated Newton or Hessian free methods, solve the Newton equation (linear system) approximately. CG method is a commonly used technique to solve the trust region subproblem approximately. In this section, we discuss inexact Newton method and its trust region variation.

3.1 Inexact Newton Method

The quadratic model $m_{k}(p)$ obtained using Taylor’s theorem is given below:

[TABLE]

Taking derivative of $m_{k}(p)$ w.r.t. $p$ and equating to zero, we get,

[TABLE]

which is Newton system and its solution gives Newton method, as given below:

[TABLE]

The computational complexity of this iteration is $O(n^{2}l+n^{3})$ which is very expensive. This iteration involves the calculation and inversion of a large Hessian matrix which is not only very expensive to calculate but expensive to store also. CG method approximately solves the subproblem (8) without forming the Hessian matrix, which solves the issues related to large computational complexity and need to store the large Hessian matrix. Each iteration runs for a given number of CG iterations or until the residual condition is satisfied, as given below:

[TABLE]

where $r_{k}=\nabla^{2}F\left(w_{k}\right)p+\nabla F\left(w_{k}\right)$ and $\eta_{k}^{{}^{\prime}}$ is a small positive value, known as forcing term Nocedal and Wright (1999).

3.2 Trust Region Inexact Newton Method

Trust region is a region in which the approximate quadratic model of the given function gives correct approximation for that function. In trust region methods, we don’t need to calculate the step size (also called learning rate) directly but they indirectly adjust the step size as per the trust region radius. Trust region method solves the following subproblem to get the step direction $p_{k}$:

[TABLE]

where $m_{k}(p)$ is a quadratic model of $F\left(w_{k}+p\right)-F\left(w_{k}\right)$, as given in (7) and $\triangle_{k}$ is the trust region radius. This subproblem can be solved similar to Newton-CG, except that now we need to take care of the extra constraint of $p$. TRON (trust region Newton method) Lin et al. (2008) is one of the most famous and widely used such method, which is used in LIBLINEAR Fan et al. (2008) to solve l2-regularized logistic regression and l2-SVM problems. Hsia et al. (2017) extends TRON by proposing better trust region radius. Hsia et al. (2018) further extends TRON using PCG subproblem solver which uses average of identity matrix and diagonal matrix as preconditioner, to solve the trust region subproblem and shows that PCG could be effective to solve ill-conditioned problems.

Then the ratio of actual and predicted reductions of the model is calculated, as given below:

[TABLE]

The parameters are updated for the $(k+1)$th iteration as given below:

[TABLE]

where $\eta_{0}>0$ is a given constant. Then the trust region radius $\triangle_{k}$ is updated as per the ratio of actual reduction and predicted reduction, and a framework for updating $\triangle_{k}$ as given in Lin and Moré (1999), is given below:

[TABLE]

where $0<\eta_{1}<\eta_{2}\leq 1$ and $0<\gamma_{1}<\gamma_{2}<1<\gamma_{3}$. If $\rho_{k}\leq\eta_{1}$ then the Newton step is considered unsuccessful and the trust region radius is shrunk. On the other hand if $\rho_{k}\geq\eta_{2}$ then the step is successful and the trust region radius is enlarged. We have implemented this framework as given in the LIBLINEAR library Fan et al. (2008) and chose the following pre-defined values for the above constants: $\eta_{0}=1e-4$, $\eta_{1}=0.25$, $\eta_{2}=0.75$, $\gamma_{1}=0.25$, $\gamma_{2}=0.5$ and $\gamma_{3}=4$.

4 STRON

STRON exploits the best of both, stochastic and batch regimes, using progressive subsampling to solve the large-scale learning problems. As stochastic gradient descent (SGD) trains faster for large-scale learning problems than gradient descent (GD) due to low computations per iteration but GD is more accurate than SGD due to batch calculations, similarly STRON takes benefit of low computation during initial iterations and as it reaches the solution region it uses batch calculations to find accurate solution like TRON. The major challenge with STRON is to decide when to switch from stochastic to full-batch regime, i.e., to tune the subsampling rate.

STRON introduces stochasticity into the trust region inexact Newton method and calculates subsampled function, gradient and Hessian values to solve the trust region subproblem, as given below:

[TABLE]

where $\nabla^{2}F_{S_{k}}\left(w_{k}\right)$ and $\nabla F_{X_{k}}\left(w_{k}\right)$ are subsampled Hessian and gradient values over the subsamples $S_{k}$ and $X_{k}$, respectively, as defined below:

[TABLE]

where subsamples are increasing, i.e., $|X_{k}|<|X_{k+1}|$, $|S_{k}|<|S_{k+1}|$ and $F_{X_{k}}$ is subsampled function value used for calculating $\rho_{k}$. STRON solves (15) approximately for given number of CG iterations or until the following residual condition is satisfied:

[TABLE]

where $r_{k}=\nabla^{2}F_{S_{k}}\left(w_{k}\right)p+\nabla F_{X_{k}}\left(w_{k}\right)$.

STRON is presented by Algorithm 1. It randomly selects subsamples $S_{k}$ and $X_{k}$ for the $k$th iteration (outer iterations). $X_{k}$ and $S_{k}$ are used for calculating the gradient and Hessian values, respectively. Then it solves the trust region subproblem using CG solver (inner iterations) which uses subsampled Hessian in calculating Hessian-vector products. CG stops when residual condition, same as (17), satisfies, it reaches maximum #CG iterations or it reaches the trust region boundary. The ratio of reduction in actual and predicted reduction is calculated similar to (12) but using subsampled function, and is used for updating the parameters as given in (13). Then trust region radius $\triangle_{k}$ is updated as per $\rho_{k}$ as given in (14) and these steps are repeated for a given number of iterations or until convergence.

STRON uses progressive subsampling, i.e., dynamic subsampling to calculate function, gradient and Hessian values, and solves the Newton system approximately. It is effective to deal with large-scale problems since it uses subsampling and solves the subproblem approximately, without forming the Hessian matrix but using only Hessian-vector products. So it handles the complexity issues related with the Newton method.

4.1 Complexity

The complexity of trust region inexact Newton (TRON) method depends on function, gradient and CG subproblem solver. This is dominated by CG subproblem solver and is given by $O\left(nl\right)\times\#\text{CG iterations},\text{ and for sparse data, }O\left(\#nnz\right)\times\#\text{CG iterations},$ where $\#nnz$ is number of non-zeros values in the dataset. For subsampled trust region inexact Newton (STRON) method, the complexity per-iteration is given by $O\left(n|S_{k}|\right)\times\#\text{CG iterations},\text{ and for sparse data, }O\left(\#nnz_{S_{k}}\right)\times\#\text{CG iterations},$ where $\#nnz_{S_{k}}$ is number of non-zeros values in the subsample $S_{k}$. Since #CG iterations taken by TRON and STRON do not differ much so the per-iteration complexity of STRON is smaller than TRON in the initial iterations and later becomes equal to TRON due to progressive subsampling, i.e., when $|S_{k}|=N$.

4.2 Analysis

STRON method uses progressive subsampling with the assumption that eventually mini-batch size becomes equal to the whole dataset, i.e., for some value of $k\geq\bar{k}>0$, STRON becomes TRON. So we can follow the theoretical results given in TRON, which itself refers the results from Lin and Moré (1999). For $k\geq\bar{k}$, we get different convergence depending upon the value of $\eta_{k}^{{}^{\prime}}$: If $\eta_{k}^{{}^{\prime}}<1$ then STRON converges Q-linearly, if $\eta_{k}^{{}^{\prime}}\rightarrow 0$ as $k\rightarrow\infty$ then STRON has Q-superlinear convergence and when $\eta_{k}^{{}^{\prime}}\leq\kappa_{0}\|\nabla F(w_{k})\|$ for $\kappa_{0}>0$ then STRON converges at quadratic rate.

5 Experimental Results

In this section, we discuss experimental settings and results. The experiments have been conducted with the bench marked binary datasets as given in the Table 1, which are available for download from LibSVM website111https://www.csie.ntu.edu.tw/c̃jlin/libsvmtools/datasets/.

We use following methods in experimentation:

TRON Hsia et al. (2018): This is a trust region inexact Newton method without any subsampling. It uses preconditioned CG method to solve the trust region subproblem and it is present in the current version of LIBLINEAR library Fan et al. (2008).

STRON: This is the proposed stochastic trust region inexact Newton method with progressive subsampling technique for gradient and Hessian calculations. It uses CG method to solve the trust region subproblem.

STRON-PCG: This is an extension of STRON using PCG for solving the trust region subproblem, as discussed in the Subsection A.1.

STRON-SVRG: This is another extension of STRON using variance reduction for subsampled gradient calculations, as discussed in Subsection A.2.

Newton-CG Byrd et al. (2011): This is stochastic inexact Newton method which uses CG method to solve the subproblem. It uses progressive subsampling similar to STRON.

SVRG-SQN Moritz et al. (2016): This is stochastic L-BFGS method with variance reduction for gradient calculations.

SVRG-LBFGS Kolte et al. (2015): This is another stochastic L-BFGS method with variance reduction. It differs from SVRG-SQN method in approach by which Hessian information is sampled.

5.1 Experimental Setup

The datasets have been divided into 80% and 20%, for training and testing datasets, respectively to plot the convergence, and 5-fold cross-validation has been used to present results in Table 2. We have used $\lambda=1/l$ for all methods because generally it gives good results and at the same time help to reduce number of tunable hyper-parameters. #CG iterations has been set to a sufficiently large value of 25, as all the inexact Newton methods use 5-10 iterations and hardly go beyond 20 iterations. Progressive batching scheme uses initial batch size of 1% for all datasets except ijcnn which uses 10% of dataset and the subsample size is increased linearly for all datasets, except covtype where exponential rate is used. Moreover we set the rate such that subsample size equals dataset size in 5 epochs because, generally second order methods reach the solution region in 4-5 epochs and converge in 8-10 epochs. For the sake of simplicity and avoid extra sampling, we take same subsample for Hessian and gradient calculations, i.e, $S_{k}=X_{k}$. Quasi-Newton methods (SVRG-SQN and SVRG-LBFGS) use mini-batch size of 10% with stochastic backtracking line search to find the step size and same size mini-batches are taken for gradient and Hessian subsampling. Memory of $M=5$ is used in quasi-Newton methods with $L=5$ as update frequency of Hessian inverse approximation for SVRG-SQN method. All the algorithms use similar exit criterion of decrease in gradient value where we calculate gradient ($\|g_{0}\|$) at $w_{0}$ and run the algorithms until $\|g_{k}\|\leq\epsilon\|g_{0}\|$, where $\|g_{k}\|$ is gradient value at kth-iteration and $\epsilon$ is the given tolerance level. Moreover all the methods are implemented in C++222experimental results can be reproduced using the library LIBS2ML Chauhan et al. (2019a). with MATLAB interface and experiments have been executed on MacBook Air (8 GB 1600 MHz DDR3, 1.6 GHz Intel Core i5 and 256GB SSD).

5.2 Comparative Study

The experiments have been performed with strongly convex and smooth l2-regularized logistic regression problem as given below:

[TABLE]

The results have been plotted as optimality ($F(w)-F(w^{*})$) versus training time (in seconds) and accuracy versus training time for high $\epsilon$(=$10^{-10}$)-accuracy solutions, as given in the Figs. 1, 2 and 3. As it is clear from the results, STRON converges faster than all other methods and shows improvement against TRON on accuracy vs. time plots. Moreover, quasi-Newton methods converges slower than inexact Newton methods as already established in the literature Lin et al. (2008). As per the intuitions, STRON takes initial advantage over TRON due to subsampling and as it reaches the solution region the progressive batching scheme reaches the full batching scheme and converges with same rate as TRON. That’s why, in most of the figures, we can observe STRON and TRON converging in parallel lines. Moreover, we observe a horizontal line for accuracy vs. time plot with covtype dataset because all methods give 100% accuracy.

Generally, the models are trained for low $\epsilon$(=$10^{-02}$)-accuracy solutions. So we present the results for such a case using Table 2, which reports results using 5-fold cross validation. As it is clear from the table, STRON either outperforms other solvers or shows results pretty close to the best method.

5.3 Results with SVM

We extend STRON to solve l2-SVM problem which is a non smooth problem, as given below:

[TABLE]

The results are reported in the Fig. 4 with news20 and real-sim datasets. As it is clear from the figure, STRON shows results similar to logistic regression problem and outperforms all other methods.

6 Conclusion

We proposed a novel stochastic trust region inexact Newton method, called as STRON, to solve the large-scale learning problems. The proposed method used progressive batching scheme to deal with noisy approximations of gradient and Hessian, and enjoyed the benefits of both stochastic and full batch regimes. STRON has been extended to use preconditioned CG as trust region subproblem solver and to use variance reduction for noisy gradient calculations. Our empirical results proved the efficacy of STRON against the state-of-art techniques with bench marked datasets.

Acknowledgements.

First author is thankful to Ministry of Human Resource Development, Government of INDIA, to provide fellowship (University Grants Commission - Senior Research Fellowship) to pursue his PhD. We are also thankful to the anonymous reviewers for their constructive comments to improve the quality of our manuscript.

Conflict of interest

The authors declare that they have no conflict of interest.

Appendix A Extensions

In this section, we discuss extensions of the proposed method with PCG for solving the trust region subproblem, and with variance reduction technique.

A.1 PCG Subproblem Solver

Number of iterations required by CG method to solve the subproblem depend on the condition number of the Hessian matrix. So for ill-conditioned problems CG method converges slowly. To avoid such situations, generally a non-singular matrix $M$, called preconditioner, is used as follow. For the linear system $\nabla^{2}F(w)p=-\nabla F(w)$, we solve following system:

[TABLE]

Generally, $M^{-1}=LL^{T}$ is taken to ensure the symmetry and positive definiteness of $M^{-1}\nabla^{2}F(w)$. PCG can be useful for solving the ill-conditioned problems but it involves extra computational overhead. We follow Hsia et al. (2018) to use PCG as a weighted average of identity matrix and diagonal matrix of Hessian, as given below:

[TABLE]

where $H$ is a Hessian matrix and $0\leq\alpha\leq 1$. For $\alpha=0$, there is no preconditioning and for $\alpha=1$ it is a diagonal preconditioner. In the experiments, we have taken $\alpha=0.01$ for TRON and STRON-PCG Hsia et al. (2018). To apply PCG to trust region subproblem, we can use Algorithm 2 without any modifications, after changing the trust region subproblem (11), as given below Steihaug (1983):

[TABLE]

where $\hat{p}=L^{T}p$. STRON using PCG as a trust region subproblem solver is denoted by STRON-PCG and the results are reported in Fig. 5. It compares TRON, STRON and STRON-PCG on news20 and rcv1 datasets. As it is clear from the figure, both STRON and STRON-PCG outperform TRON.

PCG trust region subproblem solver involves extra cost for calculating the preconditioner, and for TRON the overhead due to preconditioner is given by

[TABLE]

And for STRON-PCG, preconditioner involves extra cost as given below:

[TABLE]

A.2 Stochastic Variance Reduced Trust Region Inexact Newton Method

Recently, researchers have proposed stochastic variants of second order methods with variance reduction. So it is an interesting question to know that how will variance reduction work with stochastic trust region inexact Newton methods, as this is not studied yet. Our empirical results prove that variance reduction does not work in this case, even after using progressive subsampling for Hessian calculation.

To improve the quality of search direction, we have used SVRG as variance reduction technique for gradient calculations, as given below:

[TABLE]

where $\bar{w}$ is parameter value at the start of outer iteration. STRON-SVRG uses variance reduction for gradient calculations and progressive batching for Hessian calculations, as given in the Algorithm 3.

The experimental results are presented in Fig. 6 with news20 and rcv1 datasets. As it is clear from the figures, STRON-SVRG lags behind STRON and TRON, i.e., variance reduction in STRON-SVRG is not sufficient to beat the progressive batching in gradient calculations of STRON. This is because both, STRON-SVRG and STRON, are stochastic/subsampled variants of TRON and to compensate for the noisy gradient calculations, former uses well-known variance reduction strategy but later uses progressive subsampling strategy. STRON is able to beat TRON but STRON-SVRG fails and lags behind both.

Appendix B More Results

In this appendix, we provide more results and study the effect of regularization coefficient on the methods.

B.1 More Experiments

Fig. 7 presents more results on gisette and webspam datasets on l2-SVM and l2-regularized logistic regression, respectively. We observe results similar to Figs. 1–4, which show that STRON outperforms all other techniques.

B.2 Effect of Regularization Coefficient

Here we study the effect of the value of regularization coefficient ($\lambda$) on the convergence and accuracy of STRON and TRON methods. Fig. 8 presents the results with a series of values for $\lambda~{}=~{}\{1/l,1e-1,1e-3,1e-5,1e-7\}$ using news20 dataset with SVM. From the figure, it is clear that both the methods are affected by the choice of value $\lambda$. For larger values of $\lambda~{}=~{}\{1e-1,1e-3\}$, both the methods converge to the less accurate solution, as depicted in terms of optimality gap and accuracy plots. On the other hand, for values of $\lambda<1e-3$, there do not seem be much difference in the accuracy of both the methods. But, in terms of optimality gap, it is clear that smaller the value of $\lambda$, better is the solution, although there seem to be no difference on the corresponding accuracy plot. Moreover, for all the values of $\lambda$, STRON clearly outperforms TRON method.

Generally, machine learning problems involve tuning a lot of hyper-parameters which are quite difficult to tune. So to reduce number of parameters to tune, we set $\lambda~{}=1/l$, which works well in practice, as it sets the value relative to the size of the dataset and gives good results, as is clear from the figure.