{'fraction_non_alphanumeric': 0.04107662622232397, 'fraction_numerical': 0.017987376132387087, 'mean_word_length': 4.086161011310712, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Large Scale Empirical Risk Minimization via Truncated Adaptive Newton Method

We consider large scale empirical risk minimization (ERM) problems, where both the problem dimension and variable size is large. In these cases, most second order methods are infeasible due to the high cost in both computing the Hessian over all samples and computing its inverse in high dimensions. In this paper, we propose a novel adaptive sample size second-order method, which reduces the cost of computing the Hessian by solving a sequence of ERM problems corresponding to a subset of samples and lowers the cost of computing the Hessian inverse using a truncated eigenvalue decomposition. We show that while we geometrically increase the size of the training set at each stage, a single iteration of the truncated Newton method is sufficient to solve the new ERM within its statistical accuracy. Moreover, for a large number of samples we are allowed to double the size of the training set at each stage, and the proposed method subsequently reaches the statistical accuracy of the full training set approximately after two effective passes. In addition to this theoretical result, we show empirically on a number of well known data sets that the proposed truncated adaptive sample size algorithm outperforms stochastic alternatives for solving ERM problems.


Introduction
The recent advances in large scale machine learning focus largely on solving of the expected risk minimization problem, in which a set of model parameters of dimension p are found that minimize an expected loss function. More typically, the expected loss function is taken with respect to an unknown probability distribution. Therefore, the expected risk is estimated with a statistical average over N samples, where N is very large. The minimization over the statistical average is called the empirical risk minimization (ERM) problem. Indeed the computational complexity of solving ERM problems depends on the size of the dataset N if we use deterministic methods that operate on the full dataset at each iteration. Stochastic optimization has long been used to reduce this cost. Among the stochastic optimization methods are stochastic gradient descent algorithm (Robbins and Monro, 1951;Bottou, 2010), Nesterov-based methods (Nesterov et al., 2007;Beck and Teboulle, 2009), variance reduction methods (Johnson and Zhang, 2013;Nguyen et al., 2017), stochastic average gradient algorithms (Roux et al., 2012;Defazio et al., 2014a), stochastic majorization-minimization algorithms (Defazio et al., 2014b;Mairal, 2013), hybrid methods (Konecnỳ and Richtárik, 2013), and dual coordinate methods Zhang, 2013, 2016). Although these stochastic firstorder methods succeed in reducing the computational complexity of deterministic methods, they suffer from slow convergence in ill conditioned problems. This drawback inspires the development of stochastic second order methods, which improve upon the performance of first-order methods by using a curvature correction. These include subsampled Hessian (Erdogdu and Montanari, 2015;Pilanci and Wainwright, 2015;Roosta-Khorasani and Mahoney, 2016a,b), incremental Hessian updates (Gürbüzbalaban et al., 2015), stochastic dual Newton ascent (Qu et al., 2016), and stochastic quasi-Newton methods (Schraudolph et al., 2007;Ribeiro, 2014, 2015;Lucchi et al., 2015;Moritz et al., 2016;Byrd et al., 2016;Mokhtari et al., 2017). However, many of these cannot improve upon the asymptotic convergence rate of linearly convergent first-order methods. Those that do improve rates do so at great computational cost due to the computation of the full gradient and the Hessian inverse, both computationally prohibitive when N and p are large, as is often true in real datasets.
Each of these methods solve the ERM problem using the full sample set. However, the samples in ERM problems are drawn from a common distribution, so a smaller ERM problem using a smaller subset of samples should have a solution close to that of the full problem. Solving a sequence of smaller ERM problems using a subset of samples thus may reduce computational complexity. The work in (Mokhtari et al., 2016), for instance, reduces complexity of Newton's method by adaptively increasing sample size, but remains impractical for high dimensional problems due to the inverse computation. However, a key insight in ERM problems is that only the empirical average is optimized, and not the true expected loss function. Indeed this error between the optimal empirical loss and expected loss, called the statistical accuracy, can be made arbitrarily small by increasing the number of samples N . More importantly, however, is that the unavoidable presence of such an error in fact grants us latitude in the accuracy of our optimization methods. Given that we can only optimize the expected loss function up to a certain accuracy, we can employ further approximation techniques with negligible loss, so long as the error induced by these approximations is small relative to that of the stochastic approximation. It is therefore possible to significantly reduce the computational cost of our optimization methods without negatively impacting the overall accuracy, thus making second order methods a viable option.
In this paper, we propose a novel adaptive sample size second-order method, which reduces the cost of computing the Hessian by solving a sequence of ERM problems corresponding to a subset of samples and lowers the cost of computing the inverse of the Hessian by using a truncated eigenvalue decomposition. In the presented scheme, we increase the size of the training set geometrically, while solving subproblems up to their statical accuracy. We show that we can increase the sample size in such a manner that a single iteration of a truncated Newton method is sufficient to solve the increased sample size ERM problem to within its statistical accuracy. This is achieved by using the quadratic convergence region of the new ERM problem. While the proposed method does not converge at a purely quadratic rate due to truncation of the Newton step, the additional linear term incurred by the approximation can be made negligible with respect to the statistical accuracy of the data set. The resulting k-Truncated Adaptive Newton (k-TAN) method uses a rank-k approximation of the Hessian to solve high dimensional ERM problems with large sample sizes up to the statistical accuracy at a significantly lower overall cost than existing optimization techniques. Specifically, we demonstrate that, in many cases, we can reach the statistical accuracy of the problem with a total computational cost of O((2N +(log 2 N ) log k)p 2 ) using rank k approximations of the Hessian inverse.

Problem Formulation
We consider in this paper the empirical risk minimization (ERM) problem for a convex function f (x, z), where z is a realization of a random variable Z. More specifically, we seek the optimal variable x ∈ R p that minimizes the expected loss L(x) := E Z [f (x, z)]. Define x * as the variable that minimizes the expected loss, i.e.
x * := argmin In general, the problem in (1) cannot be solved without knowing the distribution of Z. As an alternative, we traditionally consider the case that we have access to N samples of Z, labelled z 1 , z 2 , . . . , z n . Define then the functions f i (x) = f (x, z i ) for i = 1, 2, . . . , N and an associated empirical risk function L n := (1/n) n i=1 f i (x) as the statistical mean over the first n ≤ N samples. We say that function L n (x) approximates the original expected loss L(x) with statistical accuracy V n if the difference between the empirical risk function L n (x) and the expected loss L(x) is upper bounded by V n for all x with high probability (w.h.p.). The statistical accuracy V n is typically bounded by V n = O(1/ √ n) (Vapnik, 2013) or the stronger V n = O(1/n) for a set of common problems (Bartlett et al., 2006;Bottou and Bousquet, 2007).
Observe that the sampled loss function L n is of an order V n difference from the true loss function L and, consequently, any additional change of the same order has negligible effect. It is therefore common to regularize non-strongly convex loss functions L n by a term of order V n . We then seek the minimum argument of the regularized risk function R n , where c is a scalar constant. The solution x * n minimizes the regularized risk function using the first n samples, which is of order V n from the expected loss function L. It follows then that by setting n = N we find a solution x * N in (2) that solves the original problem in (1) up to the statistical accuracy of using all N samples.
The problem in (2) is strongly convex can be solved using any descent method. In particular, Newton's method uses a curvature-corrected gradient to iteratively update a variable x, and is known to converge to the optimal argument x * n at a very fast quadratic rate. To implement Newton's method, it is necessary to compute the gradient ∇R n (x) and Hessian ∇ 2 R n (x) as The variable x is updated in Newton's method as Solving (1) to the full statistical accuracy V N (i.e. solving (2) for n = N ) using Newton's method would then require the computation of individual gradients and Hessians for N functions f i for computational cost of O(N p 2 ) at each iteration. Furthermore, the computation of the Hessian inverse in (4) requires a cost of O(p 3 ), bringing at total of O(N p 2 + p 3 ) for an iteration of Newton's method using the whole dataset. For large n and p, this may become computationally infeasible. In this paper we show how this complexity can be reduced by gradually increasing the sample size n and approximating the inverse of the respective Hessian ∇ 2 R n (x).

k-Truncated Adaptive Newton (k-TAN) Method
We propose the k-Truncated Adaptive Newton (k-TAN) as a low cost alternative to solving (1) to its statistical accuracy. In the k-TAN method, at each iteration we start from a point x m within the statistical accuracy of We geometrically increase the sample size to n = αm, where α > 1, and compute x n using an approximated Newton method on the increased sample size risk function R n . More specifically, we update a decision variable x m associated with R m to a new decision variable x n associated with R n with the Newton-type update whereĤ n,k is a matrix approximating the Hessian ∇ 2 R n (x m ) and parametrized by k ∈ {1, 2, . . . , p}.
In particular we are interested in an approximationĤ n,k whose inverseĤ −1 n,k can be computed with complexity less than O(p 3 ). To define such a matrix, consider µ 1 ≥ µ 2 ≥ . . . ≥ µ p to be the eigenvalues of the Hessian of empirical risk ∇ 2 L n (x m ), with associated eigenvectors v 1 , v 2 , . . . , v p . We perform an eigenvalue decomposition of The full approximated HessianĤ n,k is subsequently defined as the rank k approximation of The inverse of the approximated HessianĤ n,k can then be computed directly using U k and Σ k aŝ Observe that setting k = p, i.e., full Hessian inverse, recovers the AdaNewton method in (Mokhtari et al., 2016). To understand how we may determine k, consider that the full Hessian computed in (3) is ∇ 2 L n (x m ) regularized by cV n I. Therefore, the eigenvalues of ∇ 2 L n (x m ) less than cV n are made negligible by the regularization, and can be left out of the approximation. We thus select the k largest eigenvalues of the Hessian which are larger than ρcV n for some truncation parameter 0 < ρ < 1.
To analyze the computational complexity of (7), observe that the inverse computation in (7) requires only the inversion of diagonal matrices, and thus the primary cost in computing the k largest eigenvalues Σ k and associated eigenvectors U k . Indeed, the truncated eigenvalue decomposition {U k , Σ k } can in general be computed with at most complexity O(kp 2 ), with recent randomized algorithms even finding {U k , Σ k } with complexity O(p 2 log k) (Halko et al., 2011). This results in a total cost of, at worst, O((log k + n)p 2 ) to perform the update in (5), thus removing a O(p 3 ) cost.
In this paper we aim to show that while we geometrically increase the size of the training set, a single iteration of the truncated Newton method in (5) is sufficient to solve the new risk function within its statistical accuracy. To state this result we first assume the following assumptions hold.

Assumption 1
The loss functions f (x, z) are convex with respect to x for all values of z. Moreover, their gradients ∇f (x, z) are Lipschitz continuous with constant M . Assumption 2 The loss functions f (x, z) are self-concordant with respect to x for all z. Assumption 3 The difference between the gradients of the empirical loss L n and the statistical average loss L is bounded by V 1/2 n for all x with high probability, Based on Assumption 1, we obtain that the regularized empirical risk gradients ∇R n are Lipschitz continuous with constant M + cV n . Assumption 2 states the loss functions are additionally self concordant which is a customary assumption in the analysis of second-order methods. It also follows that the functions R n are therefore self concordant. Assumption 3 bounds the difference between gradients of the expected loss and the empirical loss with n samples by V 1/2 n . This is a reasonable bound for the convergence of gradients to their statistical averages using the law of large numbers.
We are interested in establishing the result that, as we increase n at each step, the k-TAN method stays in the quadratic region of the the associated risk function. More explicitly, we wish to show the sample size can be increased from m to n = αm such that x m is in the quadratic region of R n . Moreover, if x m is indeed in the quadratic region of R n , then we demonstrate that a single step of k-TAN as in (5) produces a point x n that is within the statistical accuracy V n of the risk R n .
Theorem 1 Consider the k-TAN method defined in (5)-(7) and suppose that the constant k for low rank factorization is defined as k = min{k | µ k+1 ≤ ρcV n } where ρ is a free parameter chosen from the interval (0, 1]. Further consider the variable x m as a V m -optimal solution of the risk R m , i.e., a solution such that Let n = αm > m and suppose Assumptions (1)-(3) hold. If the sample size n is chosen such that the following conditions are satisfied, where K = 3 + (2 + c x * 2 /2)(1 − 1/α), then the variable x n computed from (5) has the suboptimality of V n with high probability, i.e., The result in Theorem 1 establishes the required conditions to guarantee that the iterates x n always stay within the statistical accuracy of the risk R n . The expression in (9) provides a condition on growth rate α = n/m to ensures that iterate x m , which is a V m -suboptimal solution for R m , is within the quadratic convergence neighborhood of Newton's method for R n . The second condition in (10) ensures that a single iteration of k-TAN is sufficient for the updated variable x n to be within the statistical accuracy of R n . Note that the first term in the left hand side of (10) is quadratic with respect to V m and comes from the quadratic convergence of Newton's method, while the second and third terms of respective orders V 3/2 m and V m are the outcome of Hessian approximation. Indeed, these terms depend on ρ, which is the upper bound on ratio of the discarded eigenvalues µ k+1 , . . . , µ p to the regularization cV n . The truncation must be enough such that ρ is sufficiently small to make (10) hold. It is worth mentioning, as a sanity check, if we set ρ = 0 then we will keep all the eigenvalues and recover the update of Newton's method which makes the non-quadratic terms in (10) zero.
The conditions in Theorem 1 are cumbersome but can be simplified if we focus on large m and assume that the inequality V m ≤ αV n holds for n = αm. Then, (9) and (10) respectively. Observe that first condition is dependent of α and the second condition depends on α and ρ. Thus, a pair (α, ρ) must be chosen that satisfies (12) for the result in Theorem 1 to hold. We point out one such pair as the parameters α = 2, c > 16(2 √ M +1) 2 and ρ = 9/(21 c x * 2 + 16+3). Consequently, when m is large we may double the sample size with each update in until n = N , after which we will have obtained a point After log 2 N iterations (roughly 2N samples processed), we solve the full risk function R N to within the statistical accuracy V N . At each iteration, the truncated inverse step requires cost O(p 2 log k). Computing Hessians over 2N samples requires cos O (2N p 2 ), resulting in a total complexity of O(p 2 (2N + log 2 N log k).
In practice, these may be chosen in a backtracking manner, in which the iterate x m is updated using an estimate (α, ρ) pair. If the resulting iterate x n is not in statistical accuracy R N (x N ) − R N (x * N ) ≤ V N , the increase factor α is decreased by factor β and ρ is decreased by factor δ. Since R N (x * N ) is not known in practice, the suboptimality can be upper bounded using strong convexity as The resulting method is presented in Algorithm 1. After preliminaries and initializations in Steps 1-4, the backtracking loop starts in Step 6 with the sample size increase by rate α. After computing the gradient and Hessian in Step 7 and 8, the low rank decomposition with rate k = min{k|µ k+1 ≤ Algorithm 1 k-TAN 1: Parameters: Sample size increase constants α 0 > 1, ρ 0 < 1 and 0 < β, δ < 1 2: Input: Initial sample size n = m 0 and argument x n = x m0 with ∇R n (x n ) < ( √ 2c)V n 3: while n ≤ N do {main loop} 4: Update argument and index: x m = x n and m = n. Reset factor α = α 0 , ρ = ρ 0 .

13:
until ∇R n (x n ) < ( √ 2c)V n 14: end while ρcV n } in Step 9. The k-TAN update is then performed with (7)-(5) in Steps 9 and 10. The factors α and ρ are then decreased using the backtracking parameters and the statistical accuracy condition is checked. We stress that, while ∇R n (x n ) must be computed to check the exit condition in Step 13, the gradient for these samples must be computed in any case in the following iteration, so no additional computation is added by this step.

Convergence Analysis
We study the convergence properties of the k-TAN method and in particular prove the result in Theorem 1. Namely, we show that the update described in (5) produces a variable x n that is within the statistical accuracy V n . This, in turn, implies all future updates will be within the statistical accuracy of their respective increased sample functions R n until the full set of N samples is used, at which point we will have reached a point x N that is within the statistical accuracy V N of problem (1).

Preliminaries
Before proceeding with the analysis of k-TAN, we first present two propositions that relate current iterate x m to the suboptimality and quadratic convergence region to the increased sample size risk R n . Define S n (x) to be the n-suboptimality of point x with respect to R n , i.e.
where x * n is the point that minimizes R n . We establish in the following proposition a bound on the n-suboptimality of x m from the difference in sample sizes m and n and their associated statistical accuracies. The proof can be found in (Mokhtari et al., 2016, Proposition 3).
Proposition 2 Consider a point x m that minimizes the risk function R m to within its statistical accuracy V m , i.e. S m (x m ) ≤ V m . If the sample size is increased from m to n = αm and V m = αV n , the sub-optimality S n (x m ) = R n (x m ) − R n (x * n ) is upper bounded as Proposition 2 demonstrates a bound on the n-suboptimality S n (x m ) of a point x m whose msuboptimality is within statistical accuracy V m . It is also necessary to establish conditions on increase rate α such the x m is also in the quadratic convergence region of R n . Traditional analysis of Newton's method characterizes quadratic convergence in terms of the Newton decrement λ n (x) := ∇ 2 R n (x) −1/2 ∇R n (x) . The iterate x is said to be in the quadratic convergence region of R n when λ n (x) < 1/4-see (Boyd and Vandenberghe, 2004, Chapter 9.6.4). The conditions for current iterate x m to be within this region are presented in the following proposition. The proof can be found in (Mokhtari et al., 2016, Proposition 5).
Proposition 3 Define x m as an V m optimal solution of the risk R m , i.e., 1/2 as the Newton decrement of variable x associated with the risk R n . If Assumption 1-3 hold, then Newton's method at point x m is in the quadratic convergence phase for the objective function R n , i.e., λ n (w m ) < 1/4, if we have

Analysis of k-TAN
To analyze the k-TAN method, it is necessary to study the error incurred from approximating the Hessian inverse in (6) with rank k. Because we are only interested in solving each risk function R n to within its statistical accuracy V n , however, some approximation error can be afforded. In the following Lemma, we characterize the error between an approximate and exact Newton steps using the chosen rank k of the approximation and the associated eigenvalues of the Hessian.
Lemma 4 Consider the k-TAN update in (5)-(7) for some k = {0, 1, . . . , p}. Define ǫ n := µ k+1 /(cV n ). The norm of difference in the k-TAN stepĤ −1 n,k ∇R n (x m ) and the exact Newton step H −1 n ∇R n (x m )where H n =Ĥ n,p = ∇ 2 R n (x m )-can be upper bounded as The result in Lemma 4 gives us an upper bound on the error incurred in single iteration of a rank k approximation of the Newton step versus an exact Newton step. To make ǫ n small, a sufficiently large k must be chosen such that µ k+1 is in the order of V n . The size of k will therefore depend on the distribution of the eigenvalues of particular empirical risk function. However, in practical datasets of high dimension, it is often the case that most eigenvalues of the Hessian will be close to 0, in which case k can be made very small. This trend is supported by our numerical experiments on real world data sets in Section 5 and the Appendix of this paper.
With the results of Proposition 2 and Lemma 4 in mind, we can characterize the n-suboptimality of the updated variable x n from (5). This is stated formally in the following Lemma.
Lemma 5 Consider the k-TAN update in (5)-(7). If x m is in the quadratic neighborhood of R n , i.e. λ n (x m ) < 1/4, then the n-suboptimality S n (x n ) = R n (x n ) − R(x * n ) can be upper bounded by S n (x n ) ≤ 16 (3 − ǫ n ) 4 36(1 + ǫ n ) 2 S n (x m ) 2 + 30ǫ n (1 + ǫ n )S n (x m ) 3/2 + 6ǫ 2 n S n (x m ) . With Lemma 5 we establish a bound on n-suboptimality of the x n obtained from the k-TAN update in (5). Observe that this bounded by terms proportional to the n-suboptimality of the previous point, S n (x m ). We can then establish that S n (x n ) is indeed upper bounded by the statistical accuracy V n by combing the results in (14) and (17) to obtain Theorem 1. The proofs of Theorem 1 and all supporting Lemmata are provided in the Appendix.

Experiments
We compare the performance of the k-TAN method to existing optimization methods on large scale machine learning problems of practical interest. In particular, we consider a regularized logistic loss function, with regularization parameters V n = 1/n and c = 1. The k-TAN method is compared against the second order method AdaNewton (Mokhtari et al., 2016) and two first order methods-SGD and SAGA (Defazio et al., 2014a). Here, we study the performance of these methods on the logistic regression problem for two datasets. First, the GISETTE handwritten digit classification from the NIPS 2003 feature selection challenge and, second, the well-known RCV1 dataset for classifying news stories from the Reuters database. We perform additional experiments on the ORANGE dataset for KDD Cup 2009 and the BIO dataset for KDD Cup 2004.
The GISETTE dataset includes N = 6000 samples of dimension p = 5000. We use step sizes of 0.08 for both SGD and SAGA. In both k-TAN and AdaNewton, the sample size is increased by a factor of α = 2 at each iteration (the condition ∇R n (w n ) < ( √ 2c)V n is always satisfied) starting with an initial size of m 0 = 124. For both of these methods, we initially run gradient descent on R m0 for 100 iterations so that we may begin in the statistical accuracy V m0 . For k-TAN, the truncation k is observed to be able to afford a cutoff of around 0.01p in all of our simulations.
In Figure 1, the convergence results of the four methods for GISETTE data is shown. The left plot demonstrates the sub-optimality with respect to the number of gradients, or samples, processed. In particular, k-TAN and AdaNewton compute m gradients per iterations, while SGD and SAGA compute 1 gradient per iteration. Observe that the second order methods all converge with a smaller number of total processed gradients than the first order methods, reaching after 2.5 × 10 4 samples a sub-optimality of 10 −7 . We point out that, while k-TAN only approximates the Hessian inverse, its convergence path follows that of AdaNewton exactly. Indeed, both algorithms reach the statistical accuracy of 1/N = 1.6 × 10 −4 after 15000 samples, or just over two passes over the dataset. To see the gain in terms of computation time of k-TAN over AdaNewton and other methods, we present in the right image of Figure 1 the convergence in terms of runtime. In this case, k-TAN outperforms all methods, reaching a sub-optimality of 4 × 10 −6 after 60 seconds, while AdaNewton reaches only a sub-optimality of 10 −3 after 80 seconds. Note that first order methods have lower cost per iteration than all second order methods. Thus, SAGA is able to converge to 2 × 10 −5 after 80 seconds. For a very high dimensional problem, we consider the RCV1 dataset with N = 18242 and p = 47236. We use step sizes of 0.1 for both SGD and SAGA and truncate sizes of around 0.001p for k-TAN, while keeping the parameters for the other methods the same. The results of these simulations are shown in Figure 2. In the left image, observe that, in terms of processed gradients, the second order methods again outperform the SGD and SAGA, as expected, with k-TAN again following the path of AdaNewton. Given the high dimension p, the cost of computing the inverse in AdaNewton provides a large bottleneck. The gain in terms of computation time can then be best seen in the right image of Figure 2. Observe that AdaNewton becomes entirely ineffective in this high dimension. The k-TAN method, alternatively, continues to descend at a fast rate because of the inverse truncation step. For this set k-TAN outperforms all the other methods, reaching an error of 10 −7 after 1500 seconds. Since both n and p are large, SAGA performs well on this dataset due to small cost per iteration. We perform additional numerical experiments on the ORANGE dataset used for customer relationship prediction in KDD Cup 2009. We use N = 20000 samples with dimension p = 14472. The convergence results are shown in Figure 3. Observe in the right hand plot that all second order methods, perform similarly well on this dataset. The first order methods, including SAGA, do not converge after 2000 seconds. We also note that, in this experiment, we were able to reduce the truncation size k to around 0.1% of p.
In Figure 4, we show results on the BIO dataset used for protein homology classification in KDD Cup 2004. The dimensions are N = 145751 and p = 74. In this setting, the number of samples is very large put the problem dimension is very small. Observe in Figure 4 that both k-TAN and AdaNewton greatly outperform the first order methods, due to the reduced cost in Hessian computation that comes from adaptive sample size. However, because p is small, the additional gain from the truncating in the inverse in k-TAN does not provide significant benefit relative to AdaNewton.

Discussion
We demonstrated in this paper the success of the proposed k-TAN method on solving large scale empirical risk minimization problems both theoretically and empirically. The k-TAN method reduces the total cost in solving (1) to its statistical accuracy in two ways: (i) progressively increasing the sample size to reduce the costs of computing gradients and Hessians, and (ii) using a low rank approximation of the Hessian to reduce the cost of inversion. The gain provided by k-TAN relative to existing methods is therefore most significant in large scale ERM problems with both large sample size N and dimension p. To see this, consider the alternatives previously considered • Stochastic first order methods, such as SAGA (Defazio et al., 2014a) and SVRG (Johnson and Zhang, 2013), compute a single gradient per iteration, and they have the overall complexity of O(N log(N )p) to achieve statistical accuracy of the full training set if V N = O(1/N ).
• Newton's method computes gradients and Hessians over the entire dataset and inverts a matrix of size p at each iteration, requiring a total cost of O(M (N p 2 + p 3 )), where M is number of iterations required to converge. Because Newton's method converges quadratically, M may be small, but the total cost is made large by N p 2 and p 3 .
• The AdaNewton method (Mokhtari et al., 2016) computes gradients and Hessians for a subset of the dataset and inverts a matrix of size p at each iteration, with the size of the subset increasing geometrically. By doubling the sample size every iteration, the statistical accuracy can be reached in a total of log 2 N iterations after a total of 2N passes over the dataset, for a total cost of O(2N p 2 + log 2 (N )p 3 ). While Hessian computation cost is reduced, for high dimensional problems the inversion cost of p 3 dominates and the algorithm remains costly.
The k-TAN method computes gradients and Hessians on a increasing subset of data in the same manner as AdaNewton, put reduces the inversion cost at each iteration to O(p 2 log k), resulting in a total cost of O(2N p 2 + p 2 log 2 N log k), or an effective cost of O (N p 2 ), if the size of the initial training set is large enough. For ill-conditioned problems, this method is a more feasible option as a second-order method than Newton's method or AdaNewton. This theoretical intuition is indeed supported in the empirical simulations performed on large, high dimensional datasets in this paper. Now, consider the bound for Newton decrement of the k-TAN iterate λ n (x n ) from (21). As we assume that λ n (x m ) < 1/4, we have λ n (x n ) ≤ 4 (3 − ǫ n ) 2 (1 + ǫ n )λ n (x m ) 2 + λ n (x m )ǫ n .
The bound in (31) provides us then the condition in(10) for S n (x n ) ≤ V n .