1 Introduction
The problem of optimizing a cost function expressed as the sum of a loss term over each sample in an input dataset is pervasive in machine learning. One example of a cost function of this type is the partition funtion. Partition function is a central quantity in many different learning tasks including training conditional random fields (CRFs) and loglinear models
[1], and will be of central focus in this paper. Batch methods based on the quadratic bound were recently proposed [1] for the class of problems invoving the minimization of the partition function, and performed favorably in comparison to stateoftheart techniques. This paper focuses on semistochastic extension of this recently developed optimization method. Standard learning systems based on batch methods such as BFGS and memorylimited LBFGS, steepest descent (see e.g. [2]), conjugate gradient [3] or quadratic bound majorization method [1]need to make a full pass through an entire dataset before updating the parameter vector. Even though these methods can converge quickly (sometimes in several passes through the dataset), as datasets grow in size, this learning strategy becomes increasingly inefficient. To faciliate learning on massive datasets, the community increasingly turns to stochastic methods.
Stochastic optimization methods interleave the update of parameters after only processing a small minibatch of examples (potentially as small as a single datapoint), leading to significant computational savings ([4, 5, 6]
). Due to its simplicity and low computational cost, the most popular contemporary stochastic learning technique is stochastic gradient descent (SGD)
[7, 4, 8]. SGD updates the parameter vector using the gradient of the objective function as evaluated on a single example (or, alternatively, a small minibatch of examples). This algorithm admits multiple extensions, including (i) stochastic average gradient method (SAG) that averages the most recently computed gradients for each training example [9], (ii) methods that compute the (weighted) average of all previous gradients [10, 11], (iii) averaged stochastic gradient descent method (ASGD) that computes a running average of parameters obtained by SGD [12], (iv) stochastic dual coordinate ascent, that optimizes the dual objective with respect to a single dual vector or a minibatch of dual vectors chosen uniformly at random [13, 14], (v) variance reduction techniques
[15, 16, 9, 17] (some do not require storage of gradients, c.f. [15]), (vi) majorizationminimization techniques that minimize a majoring surrogate of an objective function [18, 19] and (vii) gain adaptation techniques [20, 21].Semistochastic methods can be viewed as an interpolation between the expensive reliable updates used by full batch methods, and inexpensive noisy updates used by stochastic methods. They inherit the best of both worlds by approaching the solution more quickly when close to the optimum (like a full batch method) while simultaneously reducing the computational complexity per iteration (though less aggressively than stochastic methods). Several semistochastic extensions have been explored in previous works
[22, 23, 24]. Recently, convergence theory and sampling strategies for these methods have been explored in [25, 26] and linked to results in finite sampling theory in [27].Additionally, incorporating secondorder information (i.e. Hessian) into the optimization problem ([20, 28, 29, 30, 31, 32]) was shown to often improve the performance of traditional SGD methods which typically provide fast improvement initially, but are slow to converge near the optimum (see e.g. [25]), require stepsize tuning and are difficult to parallelize [33]. This paper focuses on semistochastic extension of a recently developed quadratic bound majorization technique [1], and we call the new algorithm semistochastic quadratic bound (SQB) method. The bound computes the update on the parameter vector using the product of the gradient of the objective function and an inverse of a secondorder term that is a descriptor of the curvature of the objective function (different than the Hessian). We discuss implementation details, in particular curvature approximation, inexact solvers, and batchsize selection strategies, which make the running time of our algorithm comparable to the gradient methods and also make the method easily parallelizable. We show global convergence of the method to a stationary point under very weak assumptions (in particular convexity is not required) and a linear convergence rate when the size of the minibatch grows sufficiently fast, following the techniques of [25]. This rate of convergence matches stateoftheart incremental techniques [15, 13, 9, 18] (furthermore it is better than in case of standard stochastic gradient methods [7, 8] which typically have sublinear convergence rate [9, 34]). Compared to other existing majorizationminimization incremental techniques [18], our approach uses much tighter bounds which, as shown in [1], can lead to faster convergence.
The paper is organized as follows: Section 2 reviews quadratic bound majorization technique. Section 3 discusses stochastic and semistochastic extensions of the bound, and presents convergence theory for the proposed methods. In particular, we discuss very general stationary convergence theory under very weak assumptions, and also present a much stronger theory, including convergence rate analysis, for logistic regression. Section 4 discusses implementation details, and Section 5 shows numerical experiments illustrating the use of the proposed methods for regularized logistic regression problems. Conclusions end the paper.
The semistochastic quadratic bound majorization technique that we develop in this paper can be broadly applied to mixture models or models that induce representations. The advantages of this technique in the batch setting for learning mixture models and other latent models, was shown in the work of [1]. In particular, quadratic bound majorization was able to find better local optima in nonconvex problems than stateofthe art methods (and in less time). While theoretical guarantees for nonconvex problems are hard to obtain, the broader convergence theory developed in this paper (finding a stationary point under weak assumptions) does carry over to the nonconvex setting.
2 Quadratic bound methods
Let
be a discrete probability space over the set of
elements, and take any loglinear density model(1) 
parametrized by a vector , where are iid inputoutput pairs, is a continuous vectorvalued function mapping and is a fixed nonnegative measure. The partition function is a scalar that ensures that is a true density, so in particular (1) integrates to :
(2) 
[1] propose a fast method to find a tight quadratic bound for , shown in the subroutine Bound Computation in Algorithm 1, which finds so that
(3) 
for any and for all .
The (regularized) maximum likelihood estimation problem is equivalent to
(4) 
where means equal up to an additive constant. The bound (3) suggests the iterative minimization scheme
(5) 
where and are computed using Algorithm 1, is the regularization term and is the step size at iteration .
In this paper, we consider applying the bound to randomly selected batches of data; any such selection we denote or .
Input Parameters and for , 
Initialize 
For each 
Subroutine Bound Computation: 
For each 
Subroutine output 
Output 
3 Stochastic and semistochastic extensions
The bounding method proposed in [1] is summarized in Algorithm 1 with at every iteration. When is large, this strategy can be expensive. In fact, computing the bound has complexity , since outer products must be summed to obtain , and each other product has complexity . When the dimension is large, considerable speedups can be gained by obtaining a factored form of , as described in Section 4.1. Nonetheless, in either strategy, the size of is a serious issue.
A natural approach is to subsample a smaller selection from the training set , so that at each iteration, we run Algorithm 1 over rather than over the full data to get . When is fixed (and smaller than ), we refer to the resulting method as a stochastic extension. If instead is allowed to grow as iterations proceed, we call this method semistochastic; these methods are analyzed in [25]. All of the numerical experiments we present focus on semistochastic methods. One can also decouple the computation of gradient and curvature approximations, using different data selections (which we call and ). We show that this development is theoretically justifiable and practically very useful.
For the stochastic and semistochastic methods discussed here, the quadratic bound property (3) does not hold for , so the convergence analysis of [1] does not immediately apply. Nonetheless, it is possible to analyze the algorithm in terms of sampling strategies for .
The appeal of the stochastic modification is that when , the complexity of Algorithm 1 to compute is much lower; and then we can still implement a (modified) iteration (5). Intuitively, one expects that even small samples from the data can give good updates for the overall problem. This intuition is supported by the experimental results, which show that in terms of effective passes through the data, SQB is competitive with state of the art methods.
We now present the theoretical analysis of Algorithm 1. We first prove that under very weak assumption, in particular using only the Lipschitz property, but not requiring convexity of the problem, the proposed algorithm converges to a stationary point. The proof technique easily carries over to other objectives, such as the ones used in maximum latent conditional likelihood problems (for details see [1]), since it relies mainly only on the sampling method used to obtain . Then we focus on problem (4), which is convex, and strictly convex under appropriate assumptions on the data. We use the structure of (4) to prove much stronger results, and in particular analyze the rate of convergence of Algorithm 1.
3.1 General Convergence Theory
We present a general global convergence theory, that relies on the Lipschitz property of the objective and on the sampling strategy in the context of Algorithm 1. The end result we show here is that any limit point of the iterates is stationary. We begin with two simple preliminary results.
Lemma 1
If every is equally likely to appear in , then .
Proof 1
Algorithm 1 returns , where . If each has an equal chance to appear in , then
Note that the hypothesis here is very weak: there is no stipulation that the batch size be of a certain size, grow with iterations, etc. This lemma therefore applies to a wide class of randomized bound methods.
Lemma 2
Denote by
the infimum over all possible eigenvalues of
over all choices of batches ( may be ). Then satisfiesProof 2
For any vector and any realization of , we have
where depends on the data. Taking the expectation over of be above inequality gives the result.
Theorem 1
For any problem of form (4), apply iteration (5), where at each iteration are obtained by Algorithm 1 for two independently drawn batches subsets selected to satisfy the assumptions of Lemma 1. Finally, suppose also that the step sizes are square summable but not summable. Then converges to a finite value, and . Furthermore, every limit point of is a stationary point of .
Theorem 1 states the conclusions of [35, Proposition 3], and so to prove it we need only check that the hypotheses of this proposition are satisfied.
Proof 3
[35] consider algorithms of the form
In the context of iteration (5), at each iteration we have
where , and is the full gradient of the regularized problem (4). We choose
We now have the following results:

Unbiased error:
(6) where the second equality is obtained by independence of the batches and , and the last equality uses Lemma 1.

Gradient related condition:
(7) 
Bounded direction:
(8) 
Bounded second moment:
By part 1, we have
(9)
The covariance matrix of is proportional to the covariance matrix of the set of individual (datapoint based) gradient contributions, and for problems of form (4) these contributions lie in the convex hull of the data, so in particular the trace of the covariance must be finite. Taken together, these results show all hypotheses of [35, Proposition 3] are satisfied, and the result follows.
Theorem (1) applies to any stochastic and semistochastic variant of the method. Note that two independent data samples and are required to prove (6). Computational complexity motivates different strategies for selecting choose and . In particular, it is natural to use larger minibatches to estimate the gradient, and smaller minibatch sizes for the estimation of the secondorder curvature term. Algorithms of this kind have been explored in the context of stochastic Hessian methods [32]. We describe our implementation details in Section 4.
3.2 Rates of Convergence for Logistic Regression
The structure of objective (4) allows for a much stronger convergence theory. We first present a lemma characterizing strong convexity and Lipschitz constant for (4). Both of these properties are crucial to the convergence theory.
Lemma 3
The objective in (4) has a gradient that is uniformly norm bounded, and Lipschitz continuous.
Proof 4
The function has a Lipschitz continuous gradient if there exists an such that
holds for all . Any uniform bound for is a Lipschitz bound for . Define
and note . Let be the empirical density where the probability of observing is given by . The gradient of (4) is given by
(10) 
It is straightforward to check that the Hessian is given by
(11) 
where denotes the covariance matrix with respect to the empirical density function . Therefore a global bound for the Lipschitz constant is given by , which completes the proof.
Note that is strongly convex for any positive . We now present a convergence rate result, using results from [25, Theorem 2.2].
Theorem 2
4 Implementation details
In this section, we briefly discuss important implementation details as well as describe the comparator methods we use for our algorithm.
4.1 Efficient inexact solvers
The linear system we have to invert in iteration (5) has very special structure. The matrix returned by Algorithm 1 may be written as , where each column of is proportional to one of the vectors computed by the bound. When the dimensions of are large, it is not practical to compute the explicitly. Instead, to compute the update in iteration (5), we take advantage of the fact that
and use (computed with a simple modification to the bound method) to implement the action of . When ( is a minibatch size), the action of the transpose on a vector can be computed in , which is very efficient for small . The action of the regularized curvature approximation follows immediately. Therefore, it is efficient to use iterative minimization schemes, such as lsqr, conjugate gradient, or others to compute the updates. Moreover, using only a few iterations of these methods further regularizes the subproblems [36, 37].
It is interesting to note that even when , and is not invertible, it makes sense to consider inexact updates. To justify this approach, we present a range lemma. A similar lemma appears in [36] for a different quadratic approximation.
Lemma 4
For any , we have .
Proof 6
The matrix is formed by a sum of weighted outer products . We can therefore write
where , ( is the current iteration of the bound computation), and is a diagonal matrix with weights , where the quantities correspond to iterations in Algorithm (1). Since is in the range of by construction, it must also be the range of .
Lemma 4 tells us that there is always a solution to the linear system , even if is singular. In particular, a minimum norm solution can be found using the MoorePenrose pseudoinverse, or by simply applying lsqr or cg, which is useful in practice when the dimension is large. For many problems, using a small number of cg iterations both speeds up the algorithm and serves as additional regularization at the earlier iterations, since the (highly variable) initially small problems are not fully solved.
4.2 Minibatches selection scheme
In our experiments, we use a simple linear interpolation scheme to grow the batch sizes for both the gradient and curvature term approximations. In particular, each batch size (as a function of iteration ) is given by
where represents the cap on the maximum allowed size, is the initial batch size, and gives the rate of increase. In order to specify the selections chosen, we will simply give values for each of . For all experiments, the cap on the gradient computation was the full training set, the cap for the curvature term was taken to be 200, initial and were both set to . At each iteration of SQB, the parameter vector is updated as follows:
where ( is the step size; we use constant step size for SQB in our experiments). Notice that is the solution to the linear system and can be efficiently computed using the lsqr solver (or any other iterative solver). For all experiments, we ran a small number () iterations of lsqr, where was chosen from the set , before updating the parameter vector (this technique may be viewed as performing conjugate gradient on the bound), and chose with the best performance (the fastest and most stable convergence).
4.3 Step size
One of the most significant disadvantages of standard stochastic gradient methods [7, 8] is the choice of the step size. Stochastic gradient algorithms can achieve dramatic convergence rate if the step size is badly tuned [38, 9]. An advantage of computing updates using approximated curvature terms is that the inversion also establishes a scale for the problem, and requires minimal tuning. This is well known (phenomenologically) in inverse problems. In all experiments below, we used a constant step size; for well conditioned examples we used step size of , and otherwise .
4.4 Comparator methods
We compared SQB method with the variety of competitive stateoftheart methods which we list below:

LBFGS: limitedmemory BFGS method (quasiNewton method) tuned for loglinear models which uses both first and secondorder information about the objective function (for LBFGS this is gradient and approximation to the Hessian); we use the competitive implementation obtained from http://www.di.ens.fr/ mschmidt/Software/minFunc.html

SGD: stochastic gradient descent method with constant step size; we use the competitive implementation obtained from http://www.di.ens.fr/ mschmidt/Software/SAG.html which is analogous to L. Bottou implementation but with prespecified step size

ASGD: averaged stochastic gradient descent method with constant step size; we use the competitive implementation obtained from http://www.di.ens.fr/ mschmidt/Software/SAG.html which is analogous to L. Bottou implementation but with prespecified step size

SAG: stochastic average gradient method using the estimate of Lipschitz constant at iteration set constant to the global Lipschitz constant; we use the competitive implementation of [9] obtained from http://www.di.ens.fr/ mschmidt/Software/SAG.html

SAGls: stochastic average gradient method with line search, we use the competitive implementation of [9] obtained from http://www.di.ens.fr/ mschmidt/Software/SAG.html; the algorithm adaptively estimates Lipschitz constant
with respect to the logistic loss function using linesearch
Since our method uses the constant step size we chose to use the same scheme for the competitor methods like SGD, ASGD and SAG. For those methods we tuned the step size to achieve the best performance (the fastest and most stable convergence). Remaining comparators (LBFGS and SAGls) use linesearch.
5 Experiments
We performed experiments with regularized logistic regression on binary classification task with regularization parameter . We report the results for six datasets. The first three are sparse: rcv1 (, ; SQB parameters: , , ), adult (, ; SQB parameters: , , ) and sido (, ; SQB parameters: , , ). The remaining datasets are dense: covtype (, ; SQB parameters: , , ), protein (, ; SQB parameters: , , ) and quantum (, ; SQB parameters: , , ). Each dataset was split to training and testing datasets such that of the original datasets was used for training and the remaining part for testing. Only sido and protein were split in half to training and testing datasets due to large disproportion of the number of datapoints belonging to each class. The experimental results we obtained are shown in Figure 1. We report the training and testing costs as well as the testing error as a function of the number of effective passes through the data and thus the results do not rely on the implementation details. We would like to emphasize however that under current implementation the average running time for the bound method across the datasets is comparable to that of the competitor methods. All codes are released and are publicly available at www.columbia.edu/ aec2163/NonFlash/Papers/Papers.html.
6 Conclusions
We have presented a new semistochastic quadratic bound (SQB) method, together with convergence theory and several numerical examples. The convergence theory is divided into two parts. First, we proved convergence to stationarity of the method under weak hypotheses (in particular, convexity is not required). Second, for the logistic regression problem, we provided a stronger convergence theory, including a rate of convergence analysis.
The main contribution of this paper is to apply SQB methods in a semistochastic largescale setting. In particular, we developed and analyzed a flexible framework that allows samplebased approximations of the bound from [1] that are appropriate in the largescale setting, computationally efficient, and competitive with stateoftheart methods.
Future work includes developing a fully stochastic version of SQB, as well as applying it to learn mixture models and other latent models, as well as models that induce representations, in the context of deep learning.
References
 [1] T. Jebara and A. Choromanska. Majorization for CRFs and latent likelihoods. In NIPS, 2012.
 [2] Y. Nesterov. Introductory lectures on convex optimization : a basic course. Applied optimization. Kluwer Academic Publ., Boston, Dordrecht, London, 2004.
 [3] M.R. Hestenes and E. Stiefel. Methods of conjugate gradients for solving linear systems. Journal of Research of the National Bureau of Standards, 49(6):409–436, 1952.

[4]
L. Bottou.
Online algorithms and stochastic approximations.
In David Saad, editor,
Online Learning and Neural Networks
. Cambridge University Press, Cambridge, UK, 1998.  [5] N. Littlestone. Learning quickly when irrelevant attributes abound: A new linearthreshold algorithm. Mach. Learn., 2(4):285–318, April 1988.

[6]
F. Rosenblatt.
The perceptron: A probabilistic model for information storage and organization in the brain.
Psychological Review, 65(6):386–408, 1958.  [7] H. Robbins and S. Monro. A Stochastic Approximation Method. The Annals of Mathematical Statistics, 22(3):400–407, 1951.
 [8] L. Bottou and Y. LeCun. Large scale online learning. In NIPS, 2003.
 [9] N. Le Roux, M. W. Schmidt, and F. Bach. A stochastic gradient method with an exponential convergence rate for finite training sets. In NIPS, 2012.
 [10] Y. Nesterov. Primaldual subgradient methods for convex problems. Math. Program., 120(1):221–259, 2009.
 [11] P. Tseng. An incremental gradient(projection) method with momentum term and adaptive stepsize rule. SIAM J. on Optimization, 8(2):506–531, February 1998.
 [12] B. T. Polyak and A. B. Juditsky. Acceleration of stochastic approximation by averaging. SIAM J. Control Optim., 30(4):838–855, July 1992.
 [13] S. ShalevShwartz and T. Zhang. Proximal stochastic dual coordinate ascent. CoRR, abs/1211.2717, 2012.
 [14] S. ShalevShwartz and T. Zhang. Accelerated minibatch stochastic dual coordinate ascent. CoRR, abs/1305.2581, 2013.
 [15] R. Johnson and T. Zhang. Accelerating stochastic gradient descent using predictive variance reduction. In C.J.C. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K.Q. Weinberger, editors, NIPS, pages 315–323. 2013.
 [16] C. Wang, X. Chen, A. Smola, and E. Xing. Variance reduction for stochastic gradient optimization. In NIPS, pages 181–189. 2013.
 [17] S. ShalevShwartz and T. Zhang. Stochastic dual coordinate ascent methods for regularized loss minimization. CoRR, abs/1209.1873, 2012.
 [18] J. Mairal. Optimization with firstorder surrogate functions. In ICML (3), pages 783–791, 2013.
 [19] J. Mairal. Stochastic majorizationminimization algorithms for largescale optimization. In NIPS, pages 2283–2291. 2013.
 [20] N. N. Schraudolph. Local gain adaptation in stochastic gradient descent. In ICANN, 1999.
 [21] N. N. Schraudolph. Fast curvature matrixvector products for secondorder gradient descent. Neural Computation, 14:2002, 2002.
 [22] A. Shaprio Y. Wardi. Convergence analysis of stochastic algorithms. Mathematics of Operations Research, 21(3):615–628, 1996.
 [23] Alexander Shapiro, Tito Homem de mello, and Pii S. On the rate of convergence of optimal solutions of monte carlo approximations of stochastic programs. SIAM Journal on Optimization, 11:70–86, 2000.
 [24] A. J. Kleywegt, A. Shapiro, and T. Homemde Mello. The sample average approximation method for stochastic discrete optimization. SIAM J. on Optimization, 12(2):479–502, February 2002.
 [25] M. P. Friedlander and M. Schmidt. Hybrid deterministicstochastic methods for data fitting. SIAM J. Scientific Computing, 34(3), 2012.
 [26] R. H. Byrd, G. M. Chin, J. Nocedal, and Y. Wu. Sample size selection in optimization methods for machine learning. Math. Program., 134(1):127–155, August 2012.
 [27] A. Aravkin, M. P. Friedlander, F. Herrmann, and T. van Leeuwen. Robust inversion, dimensionality reduction, and randomized sampling. Mathematical Programming, 134(1):101–125, 2012.
 [28] N. N. Schraudolph, J. Yu, and S. Günter. A stochastic quasinewton method for online convex optimization. In AISTATS, 2007.
 [29] Y. Le Cun, L. Bottou, G. B. Orr, and K.R. Müller. Efficient backprop. In Neural Networks, Tricks of the Trade, Lecture Notes in Computer Science LNCS 1524. Springer Verlag, 1998.
 [30] A. Bordes, L. Bottou, and P. Gallinari. Sgdqn: Careful quasinewton stochastic gradient descent. J. Mach. Learn. Res., 10:1737–1754, December 2009.
 [31] J. Martens. Deep learning via hessianfree optimization. In ICML, 2010.
 [32] R. H. Byrd, G. M. Chin, W. Neveitt, and J. Nocedal. On the use of stochastic hessian information in optimization methods for machine learning. SIAM Journal on Optimization, 21(3):977–995, 2011.
 [33] Q. V. Le, J. Ngiam, A. Coates, A. Lahiri, B. Prochnow, and A. Y. Ng. On optimization methods for deep learning. In ICML, 2011.
 [34] A. Agarwal, P. L. Bartlett, P. D. Ravikumar, and M. J. Wainwright. Informationtheoretic lower bounds on the oracle complexity of stochastic convex optimization. IEEE Transactions on Information Theory, (5):3235–3249.
 [35] D. P. Bertsekas and J. N. Tsitsiklis. Gradient convergence in gradient methods with errors. SIAM J. on Optimization, 10(3):627–642, July 1999.
 [36] James Martens and Ilya Sutskever. Training deep and recurrent networks with hessianfree optimization. In Grégoire Montavon, GenevièveB. Orr, and KlausRobert Müller, editors, Neural Networks: Tricks of the Trade, volume 7700 of Lecture Notes in Computer Science, pages 479–535. Springer Berlin Heidelberg, 2012.
 [37] Curtis R. Vogel. Computational Methods for Inverse Problems. Society for Industrial and Applied Mathematics, 2002.
 [38] A. Nemirovski, A. Juditsky, G. Lan, and A. Shapiro. Robust stochastic approximation approach to stochastic programming. SIAM J. on Optimization, 19(4):1574–1609, January 2009.
Comments
There are no comments yet.