Unbiased scalable softmax optimization

TL;DR

This paper introduces the first unbiased algorithms for scalable softmax optimization that are efficient regardless of the number of classes or data points, significantly improving performance on large-scale tasks.

Contribution

It presents novel unbiased algorithms for softmax optimization with constant work per iteration, independent of class and data size, outperforming existing methods.

Findings

Unbiased algorithms outperform biased estimates on real datasets.

Work per iteration is independent of class and data size.

Methods are scalable and efficient for large neural network models.

Abstract

Recent neural network and language models rely on softmax distributions with an extremely large number of categories. Since calculating the softmax normalizing constant in this context is prohibitively expensive, there is a growing literature of efficiently computable but biased estimates of the softmax. In this paper we propose the first unbiased algorithms for maximizing the softmax likelihood whose work per iteration is independent of the number of classes and datapoints (and no extra work is required at the end of each epoch). We show that our proposed unbiased methods comprehensively outperform the state-of-the-art on seven real world datasets.

Figures

Tables

Table 1. Table 1: Datasets with a summary of their properties. Where the number of classes, dimension or number of examples has been altered, the original value is displayed in brackets.

Data set	Classes	Dimension	Examples
MNIST	10	780	60,000
Bibtex	147 (159)	1,836	4,880
Delicious	350 (983)	500	12,920
Eurlex	838 (3,993)	5,000	15,539
AmazonCat-13K	2,709 (2,919)	10,000 (203,882)	100,000 (1,186,239)
Wiki10	4,021 (30,938)	10,000 (101,938)	14,146
WikiSmall	18,207 (28,955)	10,000 (2,085,164)	90,737 (342,664)

Table 2. Table 2: Tuned initial learning rates for each algorithm on each dataset. The learning rate in 10 0 , ± 1 , ± 2 , ± 3 / N superscript 10 0 plus-or-minus 1 plus-or-minus 2 plus-or-minus 3 𝑁 10^{0,\pm 1,\pm 2,\pm 3}/N with the lowest log-loss after 50 epochs using only 10% of the data is displayed. Vanilla SGD applied to AmazonCat, Wiki10 and WikiSmall suffered from overflow with a learning rate of 10 − 3 / N superscript 10 3 𝑁 10^{-3}/N , but was stable with smaller learning rates (the largest learning rate for which it was stable is displayed).

Data set	OVE	NCE	IS	Vanilla	Umax	Implicit
MNIST	${10}^1$	${10}^1$	${10}^1$	${10}^{-2}$	${10}^1$	${10}^{-1}$
Bibtex	${10}^2$	${10}^2$	${10}^2$	${10}^{-2}$	${10}^{-1}$	${10}^1$
Delicious	${10}^1$	${10}^3$	${10}^3$	${10}^{-3}$	${10}^{-2}$	${10}^{-2}$
Eurlex	${10}^{-1}$	${10}^2$	${10}^2$	${10}^{-3}$	${10}^{-1}$	${10}^1$
AmazonCat	${10}^1$	${10}^3$	${10}^3$	${10}^{-5}$	${10}^{-2}$	${10}^{-3}$
Wiki10	${10}^{-2}$	${10}^3$	${10}^2$	${10}^{-4}$	${10}^{-2}$	${10}^0$
WikiSmall	${10}^3$	${10}^3$	${10}^3$	${10}^{-4}$	${10}^{-3}$	${10}^{-3}$

Table 3. Table 3: Relative log-loss after 50 epochs. The values for each dataset are normalized by dividing by the corresponding Implicit log-loss. The algorithm with the lowest log-loss for each dataset is in bold.

Data set	OVE	NCE	IS	Vanilla	U-max	Implicit
MNIST	5.25	5.55	5.26	1.31	1.40	1.00
Bibtex	12.65	12.65	12.48	6.61	4.25	1.00
Delicious	1.77	1.78	1.76	1.16	1.03	1.00
Eurlex	4.65	4.59	4.58	2.58	1.50	1.00
AmazonCat	2.01	2.03	2.00	1.39	0.93	1.00
Wiki10	3.68	3.72	3.64	3.13	1.24	1.00
WikiSmall	1.33	1.33	1.33	1.13	1.01	1.00
Average	4.48	4.52	4.44	2.47	1.62	1.00

Table 4. Table 4: Time in seconds taken to run 50 epochs. OVE/NCE/IS/Vanilla/U-max with n = 1 , m = 5 formulae-sequence 𝑛 1 𝑚 5 n=1,m=5 all have the same runtime. Implicit SGD with n = 1 , m = 1 formulae-sequence 𝑛 1 𝑚 1 n=1,m=1 is faster per iteration. The final column displays the ration of OVE/…/U-max to Implicit SGD for each dataset.

Data set	Implicit SGD	OVE/NCE/IS/Vanilla/U-max	Ratio
MNIST	1283	2494	1.94
Bibtex	144	197	1.37
Delicious	287	325	1.13
Eurlex	427	903	2.12
AmazonCat	24392	42816	1.76
Wiki10	783	1223	1.56
WikiSmall	6407	8470	1.32
Average	-	-	1.60

Equations

$$p(y=\ell |x;W)=\frac{e^{x^{\top }{w}_{\ell }}}{{\sum }_{k=1}^K{e}^{x^{\top }{w}_k}},$$

$$L(W)$$

$$-\log (a)$$

$$=\underset{u}{\max }\{-u-\exp (-u)a+1\}$$

$$f(u,W)=\sum _{i=1}^N\sum _{k\ne y_i}\frac{u_i+e^{-u_i}}{K-1}+e^{x_i^{\top }(w_k-w_{y_i})-u_i}+\frac{\mu }{2}\Vert W{\Vert }_2^2$$

$$f_{ik}(u,W)=N\left(u_i+e^{-u_i}+(K-1)e^{x_i^{\top }(w_k-w_{y_i})-u_i}\right)+\frac{\mu }{2}({\beta }_{y_i}\Vert w_{y_i}{\Vert }_2^2+{\beta }_k\Vert w_k{\Vert }_2^2),$$

$${\nabla }_{w_k}{f}_{ik}$$

$${\nabla }_{w_{y_i}}{f}_{ik}$$

$${\nabla }_{w_j}{f}_{ik}$$

$${\nabla }_{u_i}{f}_{ik}$$

$$L(W)=-\underset{\bar{u}}{\min }\left\{\frac{1}{N}\sum _{i=1}^N\frac{1}{K-1}\sum _{k\ne y_i}{\bar{f}}_{ik}(\bar{u},W)\right\}+N$$

$$\bar{f}_{ik}(\bar{u},W)=N\big{(}\bar{u}_{i}-x_{i}^{\top}w_{y_{i}}+e^{x_{i}^{\top}w_{y_{i}}-\bar{u}_{i}}+(K-1)e^{x_{i}^{\top}w_{k}-\bar{u}_{i}}\big{)}+\frac{\mu}{2}(\beta_{y_{i}}\|w_{y_{i}}\|_{2}^{2}+\beta_{k}\|w_{k}\|_{2}^{2})$$

$${\theta }^{(t+1)}={\theta }^{(t)}-{\eta }_t\nabla f({\theta }^{(t+1)},{\xi }_t),$$

$$\Vert \nabla f({\theta }^{(t+1)},{\xi }_t){\Vert }_2\le \Vert \nabla f({\theta }^{(t)},{\xi }_t){\Vert }_2-m\Vert {\theta }^{(t+1)}-{\theta }^{(t)}{\Vert }_2$$

$$f_{ik}(v,W)=$$

$${\nabla }_{w_k}{f}_{ik}$$

$${\nabla }_{w_{y_i}}{f}_{ik}$$

$${\nabla }_{w_j}{f}_{ik}$$

$${\nabla }_{u_i}{f}_{ik}$$

$${\nabla }_{w_k}{f}_{ik}$$

$${\nabla }_{w_{y_i}}{f}_{ik}$$

$${\nabla }_{u_i}{f}_{ik}$$

$$-N\log (K)$$

$$u_i^{\ast }$$

$$\le \log (1+\sum _{k\ne y_i}{e}^{\Vert x_i{\Vert }_2(\Vert w_k{\Vert }_2+\Vert w_{y_i}{\Vert }_2)})$$

$$\le \log (1+\sum _{k\ne y_i}{e}^{2B_x{B}_w})$$

$$=\log (1+(K-1)e^{2B_x{B}_w})$$

$$f(u,W)$$

$$=\sum _{i=1}^N{a}_i^{\top }\theta +e^{-u_i}+\sum _{k\ne y_i}{e}^{b_{ik}^{\top }\theta }+\frac{\mu }{2}\Vert W{\Vert }_2^2.$$

$${\nabla }^{2}f(\theta )$$

$$=I\cdot \min \{\exp (-B_u),\mu \}$$

$$⪰0.$$

$$B_f$$

$$\underset{\Vert W{\Vert }_2^2\le B_W^2,0\le u\le B_u}{\max }\Vert \nabla f(u,W){\Vert }_2$$

$$\le \underset{\Vert W{\Vert }_2^2\le B_W^2,0\le u\le B_u}{\max }{\mathbb{E}}_{ik}\Vert \nabla f_{ik}(u,W){\Vert }_2$$

$$\le \underset{\Vert W{\Vert }_2^2\le B_W^2,0\le u\le B_u}{\max }\underset{ik}{\max }\Vert \nabla f_{ik}(u,W){\Vert }_2.$$

$$\Vert {\nabla }_{u_i}{f}_{ik}(u,W){\Vert }_2$$

Full text

1 Introduction

Under the softmax model111Also known as the multinomial logit model. the probability that a random variable $y$ takes on the label $\ell\in\{1,...,K\}$, is given by

[TABLE]

where $x\in\mathbb{R}^{D}$ is the covariate, $w_{k}\in\mathbb{R}^{D}$ is the vector of parameters for the $k$-th class, and $W=[w_{1},w_{2},...,w_{K}]\in\mathbb{R}^{D\times K}$ is the parameter matrix. Given a dataset of $N$ label-covariate pairs $\mathcal{D}=\{(y_{i},x_{i})\}_{i=1}^{N}$, the ridge-regularized maximum log-likelihood problem is given by

[TABLE]

where $\|W\|_{2}$ denotes the Frobenius norm.

The softmax is a fundamental and ubiquitous distribution, with applications in fields such as economics and biomedicine (Rust & Zahorik, 1993; Kirkwood & Sterne, 2010; Gopal & Yang, 2013) and appears as a convex surrogate for the (hard) maximum loss in discrete optimization (Maddison et al., 2016) and network flows (Shahrokhi & Matula, 1990). This paper focusses on how to maximize (2) when $N$, $K$, $D$ are all large. Large values for $N$,$K$,$D$ are increasingly common in modern applications such as natural language processing and recommendation systems, where $N$, $K$, and $D$ can each be on the order of millions or billions (Chelba et al., 2013; Partalas et al., 2015).

A natural approach to maximizing $L(W)$ with large values for $N$, $K$ and $D$ is to use Stochastic Gradient Descent (SGD), sampling a mini-batch of datapoints each iteration. However when $K$ and $D$ are large, the $O(KD)$ cost of calculating the normalizing sum $\sum_{k=1}^{K}e^{x_{i}^{\top}w_{k}}$ in the stochastic gradients can be prohibitively expensive. Several approximations that avoid calculating the normalizing sum have been proposed to address this difficulty. These include tree-structured methods (Bengio et al., 2003; Daume III et al., 2016; Grave et al., 2016; Jernite et al., 2016), sampling methods (Bengio & Senécal, 2008; Mnih & Teh, 2012; Ji et al., 2015; Joshi et al., 2017) and self-normalization (Andreas & Klein, 2015). Alternative models such as the spherical family of losses (de Brébisson & Vincent, 2015; Vincent et al., 2015) that do not require normalization have been proposed to sidestep the issue entirely (Martins & Astudillo, 2016). Krishnapuram et al. (2005) avoid calculating the sum using a maximization-majorization approach based on lower-bounding the eigenvalues of the Hessian matrix. All222The method of Krishnapuram et al. (2005) does converge to the optimal MLE, but has $O(ND)$ runtime per iteration which is not feasible for large $N$ and $D$. of these approximations are computationally tractable for large $N$, $K$ and $D$, but are unsatisfactory in that they are biased and do not converge to the optimal $W^{\ast}=\text{argmax}\,L(W)$.

Recently333This same idea has appeared multiple times in the literature. For example (Ruiz et al., 2018) use a similar idea for variational inference of the softmax. Raman et al. (2016) showed how to recast (2) as a double-sum over $N$ and $K$. This formulation is amenable to SGD that samples only one datapoint and class in each iteration, reducing the per iteration cost to $O(D)$. However, vanilla SGD applied to this formulation is unstable in that the stochastic gradients may have high variance and a high dynamic range leading to computational overflow errors. Raman et al. (2016) deal with this instability by occasionally calculating the normalizing sum for all datapoints at a cost of $O(NKD)$. Although this achieves stability, its high cost nullifies the benefit of the cheap $O(D)$ per iteration cost.

In this paper we propose two robust unbiased SGD algorithms for optimizing double-sum formulations of the softmax likelihood. The first is an implementation of Implicit SGD, a stochastic gradient method that is known to be more stable than vanilla SGD, and yet has similar convergence properties (Toulis et al., 2016). We show that the Implicit SGD updates for the double-sum formulation can be efficiently computed using a bisection method with tight initial bounds. Furthermore, we guarantee the stability of Implicit SGD by proving that the step size is asymptotically linearly bounded (unlike vanilla SGD which is exponentially bounded). The second algorithm is a new SGD method called U-max, that is guaranteed to have bounded gradients and converges to the optimal solution of (2) for all sufficiently small learning rates. This method is particularly suited to situations where calculating simultaneous inner products is cheap (for example when using GPUs).

We compare the performance of U-max and Implicit SGD to the (biased) state-of-the-art methods for maximizing the softmax likelihood which cost $O(D)$ per iteration. Both U-max and Implicit SGD outperform all other methods. Implicit SGD has the best performance with an average log-loss 4.44 times lower than the previous state-of-the-art biased methods.

In summary, our contributions in this paper are that we:

Develop an alternative softmax double-sum formulation with gradients of smaller magnitude as compared to that in Raman et al. (2016) (Section 2). 2. 2.

Derive an efficient implementation of Implicit SGD using a bisection method, analyze its runtime and bound its step size (Section 3.1). 3. 3.

Propose the U-max algorithm to stabilize the vanilla SGD updates and prove its convergence (Section 3.2). 4. 4.

Conduct experiments showing that both U-max and Implicit SGD outperform the previous state-of-the-art, with Implicit SGD having the best performance (Section 4).

2 Convex double-sum formulation

2.1 Derivation of double-sum

In order to have an SGD method that samples both datapoints and classes each iteration, we need to represent (2) as a double-sum over datapoints and classes. We begin by rewriting (2) in a more convenient form,

[TABLE]

The key to converting (3) into its double-sum representation is to express the negative logarithm using its convex conjugate444This trick is related to the bounds given in (Gopal & Yang, 2013).:

[TABLE]

where $u=-\log(-v)$ and the optimal value of $u$ is $u^{\ast}(a)=\log(a)$. Applying (4) to each of the logarithmic terms in (3) yields $L(W)=-\min_{u\geq 0}\left\{f(u,W)\right\}+N$ where

[TABLE]

is our double-sum representation that we seek to minimize. Clearly $f$ is a jointly convex function in $u$ and $W$. The variable $u_{i}$ can be thought of as an approximation to the log-normalizer, as its optimal solution is $u_{i}^{\ast}(W)=\log(1+\sum_{k\neq y_{i}}e^{x_{i}^{\top}(w_{k}-w_{y_{i}})})\geq 0$. In Appendix B we prove that the optimal $u$ and $W$ are contained in a compact convex set and that $f$ is strongly convex within this set. Thus performing projected-SGD on $f$ is guaranteed to converge to a unique optimum with a convergence rate of $O(1/T)$ where $T$ is the number of iterations (Lacoste-Julien et al., 2012).

2.2 Instability of vanilla SGD

The challenge in optimizing $f$ using SGD is that the gradients can have very large magnitudes. Observe that ${f=\mathbb{E}_{ik}[f_{ik}]}$ where $i\sim\text{unif}(\{1,...,N\})$, $k\sim\text{unif}(\{1,...,K\}-\{y_{i}\})$ and

[TABLE]

where $\beta_{j}=\frac{N}{n_{j}+(N-n_{j})/(K-1)}$ is the inverse of the probability of class $j$ being sampled either through $i$ or $k$, and $n_{j}=|\{i:y_{i}=j\}|$. The corresponding stochastic gradient is:

[TABLE]

If $u_{i}$ is at its optimal value $u_{i}^{\ast}(W)=\log(1+\sum_{k\neq y_{i}}e^{x_{i}^{\top}(w_{k}-w_{y_{i}})})$ then $e^{x_{i}^{\top}(w_{k}-w_{y_{i}})-u_{i}}\leq 1$ and the magnitude of the $N(K-1)e^{x_{i}^{\top}(w_{k}-w_{y_{i}})-u_{i}}$ terms in the gradient are bounded by $N(K-1)\|x_{i}\|_{2}$. However if $u_{i}\ll x_{i}^{\top}(w_{k}-w_{y_{i}})$, then $e^{x_{i}^{\top}(w_{k}-w_{y_{i}})-u_{i}}\ggg 1$ and the magnitude of the gradients can become extremely large.

Extremely large gradients lead to two major problems: (a) they could lead to overflow errors and cause the algorithm to crash, (b) they result in the stochastic gradient having high variance, which leads to slow convergence555The convergence rate of SGD is inversely proportional to the second moment of its gradients (Lacoste-Julien et al., 2012).. In Section 4 we show that these problems occur in practice and make vanilla SGD both an unreliable and inefficient method666The same problems arise if we approach optimizing (3) via stochastic composition optimization (Wang et al., 2016). As is shown in Appendix C, stochastic composition optimization yields near-identical expressions for the stochastic gradients in (7) and has the same stability issues..

The sampled softmax optimizers in the literature (Bengio & Senécal, 2008; Mnih & Teh, 2012; Ji et al., 2015; Joshi et al., 2017) do not have the issue of large magnitude gradients. Their gradients are bounded by $N(K-1)\|x_{i}\|_{2}$ since their approximations ensure that $u_{i}^{\ast}(W)>x_{i}^{\top}(w_{k}-w_{y_{i}})$. For example, in one-vs-each (Titsias, 2016), $u_{i}^{\ast}(W)$ is approximated by $\log(1+e^{x_{i}^{\top}(w_{k}-w_{y_{i}})})>x_{i}^{\top}(w_{k}-w_{y_{i}})$. However, since these methods only approximate $u_{i}^{\ast}(W)$, the iterates do converge to the optimal $W^{\ast}$.

The goal of this paper is to design reliable and efficient SGD algorithms for optimizing the double-sum formulation in (5). We propose two such methods: Implicit SGD (Section 3.1) and U-max (Section 3.2). But before we introduce these methods we should establish that (5) is a good choice for the double-sum formulation.

2.3 Choice of double-sum formulation

The double-sum in (5) is different to that of Raman et al. (2016). Their formulation can be derived by applying the convex conjugate substitution to (2) instead of (3). The resulting equations are

[TABLE]

where

[TABLE]

and the optimal solution for $\bar{u}_{i}$ is $\bar{u}_{i}^{\ast}(W^{\ast})=\log(\sum_{k=1}^{K}e^{x_{i}^{\top}w_{k}^{\ast}})$. The only difference between the formulations is the reparameterization $\bar{u}_{i}=u_{i}+x_{i}^{\top}w_{y_{i}}$.

Although either double-sum formulations can be used as a basis for SGD, our formulation in (5) tends to have smaller magnitude stochastic gradients, and hence faster convergence. To see this on a high level, note that typically ${x_{i}^{\top}w_{y_{i}}=\mbox{argmax}_{k}\{x_{i}^{\top}w_{k}\}}$ and so the $\bar{u}_{i}$, $x_{i}^{\top}w_{y_{i}}$ and $e^{x_{i}^{\top}w_{y_{i}}-\bar{u}_{i}}$ terms are of the greatest magnitude in (8). Although at optimality these terms should roughly cancel, this will not be the case during the early stages of optimization, leading to stochastic gradients of large magnitude. In contrast, the function $f_{ik}$ in (6) only has $x_{i}^{\top}w_{y_{i}}$ appearing as a negative exponent, and so if $x_{i}^{\top}w_{y_{i}}$ is large then the magnitude of the stochastic gradients will be small. A more rigorous version of this argument is presented in Appendix A and in Section 4 we present numerical results confirming that our double-sum formulation leads to faster convergence.

3 Stable SGD methods

3.1 Implicit SGD

One method that solves the large gradient problem is Implicit SGD777Also known to as an “incremental proximal algorithm” (Bertsekas, 2011) or “stochastic proximal iteration” (Ryu & Boyd, 2014). (Bertsekas, 2011; Ryu & Boyd, 2014; Toulis & Airoldi, 2015; Toulis et al., 2016). Implicit SGD uses the update equation

[TABLE]

where $\theta^{(t)}$ is the value of the $t^{th}$ iterate, $f$ is the function we seek to minimize and $\xi_{t}$ is a random variable controlling the stochastic gradient such that $\nabla f(\theta)=\mathbb{E}_{\xi_{t}}[\nabla f(\theta,\xi_{t})]$. The update (9) differs from vanilla SGD in that $\theta^{(t+1)}$ appears on both the left and right side of the equation, whereas in vanilla SGD it appears only on the left side. In our case $\theta=(u,W)$ and $\xi_{t}=(i_{t},k_{t})$ with $\nabla f(\theta^{(t+1)},\xi_{t})=\nabla f_{i_{t},k_{t}}(u^{(t+1)},W^{(t+1)})$.

Although Implicit SGD has similar convergence rates to vanilla SGD, it has other properties that can make it preferable over vanilla SGD. It is more robust to the learning rate (Toulis et al., 2016), which important since a good value for the learning rate is never known a priori, and is provably more stable (Ryu & Boyd, 2014, Section 5). Another property, which is of particular interest to our problem, is that it has smaller step sizes.

Proposition 1**.**

Consider applying Implicit SGD to optimizing $f(\theta)=\mathbb{E}_{\xi}[f(\theta,\xi)]$ where $f(\theta,\xi)$ is $m$-strongly convex for all $\xi$. Then

[TABLE]

and so the Implicit SGD step size is smaller than that of vanilla SGD.

Proof.

The proof is provided in Appendix D. ∎

The bound in Proposition 1 can be tightened for our particular problem. Unlike vanilla SGD whose step size magnitude is exponential in $x_{i}^{\top}(w_{k}-w_{y_{i}})-u_{i}$, as shown in (7), for Implicit SGD the step size is asymptotically linear in $x_{i}^{\top}(w_{k}-w_{y_{i}})-u_{i}$. This effectively guarantees that Implicit SGD cannot suffer from computational overflow.

Proposition 2**.**

Consider the Implicit SGD algorithm where in each iteration only one datapoint $i$ and one class $k\neq y_{i}$ is sampled. The magnitude of its step size in $W$ is ${O(x_{i}^{\top}(\frac{w_{k}}{1+\eta\mu\beta_{k}}-\frac{w_{y_{i}}}{1+\eta\mu\beta_{y_{i}}})-u_{i})}$.

Proof.

The proof is provided in Appendix E.2. ∎

The major difficulty in applying Implicit SGD is that in each iteration one has to compute a solution to (9) (Ryu & Boyd, 2014, Section 6). The tractability of this procedure is problem dependent. We show that computing a solution to (9) is indeed tractable for the problem considered in this paper. The details are laid out in full in Appendix E.

Proposition 3**.**

Consider the Implicit SGD algorithm where in each iteration $n$ datapoints and $m$ classes are sampled. The Implicit SGD update $\theta^{(t+1)}$ can be computed to within $\epsilon$ accuracy in runtime $O(n^{2}(n+m)\log(\epsilon^{-1})+nmD)$.

Proof.

The proof is provided in Appendix E.4. ∎

In Proposition 3 the $\log(\epsilon^{-1})$ factor comes from applying a first order method to solve the strongly convex Implicit SGD update equation. It may be the case that performing this optimization is more expensive than the $O(nmD)$ cost of computing the $x_{i}^{\top}w_{k}$ inner products, and so each iteration of Implicit SGD may be significantly slower than that of vanilla SGD.

Fortunately, in certain cases we can improve the runtime of solving the implicit update. If $n=1$ and we just sample one datapoint per iteration then it is possible to reduce the update to solving just a univariate strongly convex optimization problem (see Appendix E.3 for details). Furthermore, when $m=1$ and only one class is sampled per iteration then we can derive upper and lower bounds on the one-dimensional variate to be optimized over. The optimization problem can then be solved using a bisection method, with an explicit upper bound on its cost.

Proposition 4**.**

Consider the Implicit SGD algorithm with learning rate $\eta$ where in each iteration only one datapoint $i$ and one class $k\neq y_{i}$ is sampled. The Implicit SGD iterate $\theta^{(t+1)}$ can be computed to within $\epsilon$ accuracy with only two $D$-dimensional vector inner products and at most $\log_{2}(\epsilon^{-1})+\log_{2}(|x_{i}^{\top}(\frac{w_{k}}{1+\eta\mu\beta_{k}}-\frac{w_{y_{i}}}{1+\eta\mu\beta_{y_{i}}})-u_{i}|+2\eta N\|x_{i}\|_{2}^{2}+\log(2K))$ bisection method function evaluations.

Proof.

The proof is provided in Appendix E.1 and the pseudocode is presented in Algorithm 1 in the appendix. ∎

For any reasonably large dimension $D$, the cost of the two $D$-dimensional vector inner-products will outweigh the cost of the bisection, and Implicit SGD with $n=m=1$ will have roughly the same speed per iteration as vanilla SGD with $n=m=1$. This is empirically confirmed for seven real-world datasets in Section 4.1.

However, if calculating inner products is relatively cheap (for example if $D$ is small or GPUs are used), then Implicit SGD will be slower than vanilla SGD. The U-max algorithm, presented next, is stable in the same way Implicit SGD is but has the same runtime as vanilla SGD. This makes U-max an ideal choice when inner products are cheap.

3.2 U-max method

As explained in Section 2.2, vanilla SGD has large gradients when $u_{i}\ll x_{i}^{\top}(w_{k}-w_{y_{i}})$. This can only occur when $u_{i}$ is less than its optimum value for the current $W$, since $u_{i}^{\ast}(W)=\log(1+\sum_{j\neq y_{i}}e^{x_{i}^{\top}(w_{k}-w_{y_{i}})})\geq x_{i}^{\top}(w_{k}-w_{y_{i}})$. A simple remedy is to set $u_{i}=\log(1+e^{x_{i}^{\top}(w_{k}-w_{y_{i}})})$ whenever $u_{i}\ll x_{i}^{\top}(w_{k}-w_{y_{i}})$. Since $\log(1+e^{x_{i}^{\top}(w_{k}-w_{y_{i}})})>x_{i}^{\top}(w_{k}-w_{y_{i}})$ this guarantees that $u_{i}>x_{i}^{\top}(w_{k}-w_{y_{i}})$ and so the gradients will be bounded. It also brings $u_{i}$ closer888Since $u_{i}<x_{i}^{\top}(w_{k}-w_{y_{i}})<\log(1+e^{x_{i}^{\top}(w_{k}-w_{y_{i}})})<\log(1+\sum_{j\neq y_{i}}e^{x_{i}^{\top}(w_{k}-w_{y_{i}})})=u_{i}^{\ast}(W)$. to its optimal value for the current $W$ and thereby decreases the the objective $f(u,W)$.

This is exactly the mechanism behind the U-max algorithm — see Algorithm 2 in Appendix F for its pseudocode. U-max is the same as vanilla SGD except for two modifications: (a) $u_{i}$ is set equal to $\log(1+e^{x_{i}^{\top}(w_{k}-w_{y_{i}})})$ whenever $u_{i}\leq\log(1+e^{x_{i}^{\top}(w_{k}-w_{y_{i}})})-\delta$ for some threshold $\delta>0$, (b) $u_{i}$ is projected onto $[0,B_{u}]$, and $W$ onto $\{W:\|W\|_{2}\leq B_{W}\}$, where $B_{u}$ and $B_{W}$ are set so that the optimal $u^{\ast}_{i}\in[0,B_{u}]$ and the optimal $W^{\ast}$ satisfies $\|W^{\ast}\|_{2}\leq B_{W}$. See Appendix B for more details on how to set $B_{u}$ and $B_{W}$.

Proposition 5**.**

Suppose $B_{f}\geq\max_{ik}\|\nabla f_{ik}(u,W)\|_{2}$ for all $\|W\|_{2}^{2}\leq B_{W}^{2}$ and $0\leq u\leq B_{u}$. Suppose the learning rate $\eta_{t}\leq\delta^{2}/(4B_{f}^{2})$, then U-max with threshold $\delta$ converges to the optimum of (2), and the rate of convergence is at least as fast as SGD with the same learning rate.

Proof.

The proof is provided in Appendix G. ∎

U-max directly resolves the problem of extremely large gradients. Modification (a) ensures that $\delta\geq x_{i}^{\top}(w_{k}-w_{y_{i}})-u_{i}$ (otherwise $u_{i}$ would be increased to $\log(1+e^{x_{i}^{\top}(w_{k}-w_{y_{i}})})$) and so the magnitude of the U-max gradients are bounded above by $N(K-1)e^{\delta}\|x_{i}\|_{2}$.

In U-max there is a trade-off between the gradient magnitude and learning rate that is controlled by $\delta$. For Proposition 5 to apply we require that the learning rate $\eta_{t}\leq\delta^{2}/(4B_{f}^{2})$. A small $\delta$ yields small magnitude gradients, which makes convergence fast, but necessitates a small $\eta_{t}$, which makes convergence slow.

As presented above U-max only samples one datapoint and class per iteration, but it trivially generalizes to multiple datapoints and classes. If $n$ datapoints and $m$ classes are sampled, then the runtime is $O(nmD)$, due to the vector inner-product calculations. This is the same runtime as vanilla SGD as well as the state-of-the-art biased methods such as Noise Contrastive Estimation (Mnih & Teh, 2012), Importance Sampling (Bengio & Senécal, 2008) and One-Vs-Each (Titsias, 2016).

4 Experiments

Two sets of experiments were conducted to assess the performance of the proposed methods. The first compares U-max and Implicit SGD to the state-of-the-art over seven real world datasets. The second investigates the difference in performance between the two double-sum formulations discussed in Section 2.3. We begin by specifying the experimental setup and then move onto the results.

4.1 Experimental setup

Data. We used the MNIST, Bibtex, Delicious, Eurlex, AmazonCat-13K, Wiki10, and WikiSmall datasets999All of the datasets were downloaded from http://manikvarma.org/downloads/XC/XMLRepository.html, except WikiSmall which was obtained from http://lshtc.iit.demokritos.gr/., the properties of which are summarized in Table 1. Most of the datasets are multi-label and, as is standard practice (Titsias, 2016), we took the first label as being the true label and discarded the remaining labels. To make the computation more manageable, we truncated the number of features to be at most 10,000 and the training and test size to be at most 100,000. If, as a result of the dimension truncation, a datapoint had no non-zero features then it was discarded. The features of each dataset were normalized to have unit $L_{2}$ norm. All of the datasets were pre-separated into training and test sets. We only focus on the performance on the algorithms on the training set, as the goal in this paper is to investigate how best to optimize the softmax likelihood, which is given over the training set.

Algorithms. We compared our algorithms to the state-of-the-art methods for optimizing the softmax which have runtime $O(D)$ per iteration101010Raman et al. (2016) have runtime $O(NKD)$ per epoch, which is equivalent to $O(KD)$ per iteration. This is a factor of $K$ slower than the methods we compare against. In most of our experiments, the second epoch of Raman would not have even started by the time our algorithms have already nearly converged.. The competitors include Noise Contrastive Estimation (NCE) (Mnih & Teh, 2012), Importance Sampling (IS) (Bengio & Senécal, 2008) and One-Vs-Each (OVE) (Titsias, 2016). Note that these methods are all biased and will not converge to the optimal softmax MLE, but, perhaps, something close to it. For these algorithms we set $n=100,m=5$, which are standard settings111111We also experimented setting $n=1,m=5$ in these methods and there was virtually no difference in performance except the runtime was slower. . For Implicit SGD we chose to implement the version in Proposition 4 which has $n=1,m=1$ and used Brent’s method as out bisection method solver.

For U-max and vanilla SGD we set $n=1,m=5$ and for U-max the threshold parameter $\delta=1$. For both methods we also experimented with $m=1$ but obtained significantly better performance with $m=5$. The probable reason is that having a larger $m$ value decreases the variance of the gradients, making the algorithms more stable with higher learning rates and thereby improving convergence.

The ridge regularization parameter $\mu$ was set to zero and the classes were sampled uniformly for all algorithms.

Epochs, losses and runtimes. Each algorithm was run for $50$ epochs on each dataset. The learning rate was decreased by a factor of 0.9 each epoch. Both the prediction error and log-loss (2) were recorded at the end of 10 evenly spaced epochs over the 50 epochs.

The OVE, NCE, IS, Vanilla and U-max algorithms have virtually the same runtime per iteration and so their relative performance can be gauged by plotting their log-loss over the epochs. Since Implicit SGD has to solve an inner optimization problem each iteration, its runtime will be slower than that of other algorithms with $n=1,m=1$, but may be faster than algorithms with $n=1,m>1$. Thus plotting its performance over the epochs may yield an inaccurate comparison to the other algorithms with respect to runtime.

To investigate this we measured the runtime of Implicit SGD with $n=m=1$ vs vanilla SGD121212Any of OVE, NCE, IS, Vanilla and U-max could have been used since their runtimes are virtually identical. with $n=1,m=5$ for 50 epochs on each dataset. To make the runtime comparison as fair as possible, both algorithms were coded in a standard NumPy framework. The runtime of Implicit SGD is $0.65\pm 0.15$ times that of vanilla SGD (see Table 4 in Appendix H for runtimes.)131313As noted above, vanilla SGD with $m=5$ performed significantly better than with $m=1$, thus we compare to the $m=5$ runtime. The runtime of Implicit SGD was on average $1.04\pm 0.07$ times that of vanilla SGD with $n=m=1$ for both methods.. Although these results are data, implementation and hardware dependent, they strongly indicate that Implicit SGD with $n=m=1$ is faster than vanilla SGD (or any similar method) with $n=1,m=5$. Thus plotting the log-loss over the epochs gives a conservative estimate of Implicit SGD’s relative performance with respect to runtime.

Learning rate. The magnitude of the gradient differs in each algorithm, due to either under- or over-estimating the normalizing constant from (2). To set a reasonable learning rate for each algorithm on each dataset, we ran them on 10% of the training data with initial learning rates141414The learning rates are divided by $N$ to counter the stochastic gradient being proportional to $N$ and thereby make the step size independent of $N$. $\eta=10^{0,\pm 1,\pm 2,\pm 3}/N$ . The learning rate with the best performance after 50 epochs is then used when the algorithm is applied to the full dataset. The tuned learning rates are presented in Table 2. Note that vanilla SGD requires a very small learning rate, otherwise it suffered from overflow. On average the tuned vanilla SGD learning rate is 3,019 times smaller than Implicit SGD’s and 319 times smaller than U-max’s.

4.2 Results

Comparison to state-of-the-art. Plots of the performance of the algorithms on each dataset are displayed in Figure 1 with the relative performance compared to Implicit SGD given in Table 3.

The Implicit SGD method has the best performance on all datasets but one. After just one epoch its performance is better than all of the state-of-the-art biased methods are after 50. Not only does it converge faster in the first few epochs, it also converges to the optimal MLE (unlike the biased methods that prematurely plateau). On average after 50 epochs Implicit SGD’s log-loss is a factor of 4.44 times lower than that of the biased methods.

Out of the algorithms that sample more than one class per iteration, U-max’s performance is the best. It is the only algorithm to outperform Implicit SGD on a dataset (AmazonCat). Vanilla SGD’s performance is better than the previous state-of-the-art but is generally worse than U-max. The difference in performance between vanilla SGD and U-max can largely be explained by vanilla SGD requiring a smaller learning rate to avoid computational overflow.

The sensitivity of each method to the initial learning rate can be seen in Figure 2, where the results of running each method on the Eurlex dataset with learning rates $\eta=10^{0,\pm 1,\pm 2,-3}/N$ is presented. The results agree with those in Figure 1, with Implicit SGD having the best performance for most learning rate settings. This is consistent with the theoretical results proving that Implicit SGD is robust to the learning rate (Ryu & Boyd, 2014; Toulis & Airoldi, 2015). In fact, Implicit SGD’s worst performance is still better than the best performance all of the other algorithms.

For learning rates $\eta=10^{1,2}/N$ the U-max log-loss is extremely large. This can be explained by Proposition 5, which does not guarantee convergence for U-max if the learning rate is too high. Vanilla SGD only has one line plotted, corresponding to the learning rate of $10^{-3}/N$, as for any high learning rate the algorithm suffered from computational overflow. The OVE, NCE and IS methods are very robust to the learning rate, which is perhaps why they have been so popular in the past.

Comparison of double-sum formulations. Figure 3 illustrates the performance on the Eurlex dataset of U-max using the proposed double-sum in (6) compared to U-max using the double-sum of Raman et al. (2016) in (8). The proposed double-sum outperforms for all151515The learning rates $\eta=10^{1,2,3,4}/N$ are not displayed in the Figure 3 for visualization purposes. They have similar behavior as $\eta=1.0/N$. learning rates $\eta=10^{0,\pm 1,\pm 2,\pm 3,\pm 4}/N$, with its $50^{th}$-epoch log-loss being $3.08$ times lower on average. This supports the argument from Section 2.3 that SGD methods applied to the proposed double-sum have smaller magnitude gradients and converge faster. Indeed, if the log-loss of vanilla SGD, U-max and Implicit SGD in Figure 1 and Table 3 were multiplied by 3.08 they would be roughly the same as OVE, NCE and IS. Thus our proposed double-sum formulation is crucial to the success of the U-max and Implicit SGD algorithms.

5 Conclusion

In this paper we propose two unbiased robust algoritms for for optimizing the softmax likelihood: Implicit SGD and U-max. These are the first unbiased algorithms that require only $O(D)$ computation per iteration, and no additional work at the end of each epoch. Implicit SGD can be efficiently implemented and clearly out-performs the previous state-of-the-art on seven real world datasets. The result is a new method that enables optimizing the softmax for extremely large number of samples and classes.

One limitation of the Implicit SGD method is that it is relatively slow if multiple datapoints are sampled each iteration or multiple inner-products can be efficiently computed (e.g. using GPUs). U-max should be the method of choice in such a setting.

We only tested U-max and Implicit SGD on the simple softmax, but these methods can also be applied to any neural network where the final layer is the softmax. Furthermore, applying these methods to word2vec type models, which can be viewed as a softmax where both $x$ and $w$ are parameters to be fit, might be provide a significant speed-up.

Appendix A Comparison of double-sum formulations

In Section 2.3 our double-sum formulation was compared to that of Raman et al. (2016). It was noted that the formulations only differ by a reparameterization $\bar{u}_{i}=u_{i}+x_{i}^{\top}w_{y_{i}}$, and an intuitive argument was given as to why our formulation leads to smaller magnitude gradients. Here we flesh out that argument and also explore different reparameterizations.

Let us introduce the set of parameterizations $v_{i}=\log(1+\sum_{k\neq y_{i}}e^{x_{i}^{\top}(w_{k}-w_{y_{i}})})+\alpha x_{i}^{\top}w_{y_{i}}$ where $\alpha\in\mathbb{R}$. Our double-sum corresponds to $\alpha=0$ while that of Raman et al. (2016) to $\alpha=1$. The question is, what is the optimal $\alpha$? The stochastic functions with $v_{i}$ are of the form

[TABLE]

where for notational simplicity we have set the ridge-regularization parameter $\mu=0$. The stochastic gradients are

[TABLE]

Let $x_{i}^{\top}w_{y_{i}}=x_{i}^{\top}\tilde{w}_{y_{i}}+\delta_{i}$ where $\tilde{w}_{y_{i}}$ is the old value of $w_{y_{i}}$ from the previous time datapoint $i$ was sampled. Let us write $v_{i}=\tilde{u}_{i}-\gamma_{i}+\alpha(x_{i}^{\top}\tilde{w}_{y_{i}}-\epsilon_{i})$, where $\gamma_{i}$ is the error between $\log(1+\sum_{k\neq y_{i}}e^{x_{i}^{\top}(\tilde{w}_{k}-\tilde{w}_{y_{i}})})$ and its estimate $\tilde{u}_{i}$, while $\epsilon_{i}$ is the error from estimating the $x_{i}^{\top}\tilde{w}_{y_{i}}$ term. The gradients become

[TABLE]

The goal is for the variance of these stochastic gradients to be as small as possible. This may be achieved by setting $\alpha$ to decrease the effect of the noise factors $\delta_{i}$ and $\epsilon_{i}$. The noise $\epsilon_{i}$ always appears as $e^{\alpha\epsilon_{i}}$ and so it is best to have $\alpha=0$ from this perspective. The noise $\delta_{i}$ appears as $e^{\alpha\delta_{i}}$, $e^{(1-\alpha)\delta_{i}}\cdot(K-1)e^{x_{i}^{\top}(w_{k}-\tilde{w}_{y_{i}})}$ or $(1-\alpha)e^{(1-\alpha)\delta_{i}}\cdot(K-1)e^{x_{i}^{\top}(w_{k}-\tilde{w}_{y_{i}})}$, and so there is tension between setting $\alpha=0$ or $\alpha=1$.

The optimal value of $\alpha$ clearly depends on the data and algorithm. If the noise $\epsilon$ is large and it is often the case that $x_{i}^{\top}\tilde{w}_{y_{i}}>x_{i}^{\top}\tilde{w}_{k}+\log(K-1)$ then $\alpha\approx 0$ is ideal, but if not then it is best for $\alpha\approx 1$. In Section 4 we showed that for our datasets $\alpha=0$ yields better results than $\alpha=1$, although the optimal value of $\alpha$ is probably between [math] and $1$. A future line of work is to develop methods to learn the optimal $\alpha$, perhaps dynamically per datapoint.

Appendix B Proof of variable bounds and strong

convexity

We first establish that the optimal values of $u$ and $W$ are bounded. Next, we show that within these bounds the objective is strongly convex and its gradients are bounded.

Lemma 1** (Raman et al. (2016)).**

The optimal value of $W$ is bounded as $\|W^{\ast}\|_{2}^{2}\leq B_{W}^{2}$ where $B_{W}^{2}=\frac{2}{\mu}N\log(K)$.

Proof.

[TABLE]

Rearranging gives the desired result. ∎

Lemma 2**.**

The optimal value of $u_{i}$ is bounded as $u_{i}^{\ast}\leq B_{u}$ where $B_{u}=\log(1+(K-1)e^{2B_{x}B_{w}})$ and $B_{x}=\max_{i}\{\|x_{i}\|_{2}\}$

Proof.

[TABLE]

∎

Lemma 3**.**

If $\|W\|_{2}^{2}\leq B_{W}^{2}$ and $u_{i}\leq B_{u}$ then $f(u,W)$ is strongly convex with convexity constant greater than or equal to $\min\{\exp(-B_{u}),\mu\}$.

Proof.

Let us rewrite $f$ as

[TABLE]

where $\theta=(u^{\top},w_{1}^{\top},...,w_{k}^{\top})\in\mathbb{R}^{N+KD}$ with $a_{i}$ and $b_{ik}$ being appropriately defined. The Hessian of $f$ is

[TABLE]

where $e_{i}$ is the $i^{th}$ canonical basis vector, $0_{N}$ is an $N$-dimensional vector of zeros and $1_{KD}$ is a $KD$-dimensional vector of ones. It follows that

[TABLE]

∎

Lemma 4**.**

If $\|W\|_{2}^{2}\leq B_{W}^{2}$ and $u_{i}\leq B_{u}$ then the 2-norm of both the gradient of $f$ and each stochastic gradient $f_{ik}$ are bounded by

[TABLE]

Proof.

By Jensen’s inequality

[TABLE]

Using the results from Lemmas 1 and 2 and the definition of $f_{ik}$ from (6),

[TABLE]

and for $j$ indexing either the sampled class $k\neq y_{i}$ or the true label $y_{i}$,

[TABLE]

Letting

[TABLE]

we have

[TABLE]

In conclusion:

[TABLE]

∎

Appendix C Stochastic Composition Optimization

We can write the equation for $L(W)$ from (3) as (where we have set $\mu=0$ for notational simplicity),

[TABLE]

where $i\sim unif(\{1,...,N\})$, $k\sim unif(\{1,...,K\})$, $h_{i}(v)\in\mathbb{R}$, $g_{k}(W)\in\mathbb{R}^{N}$ and

[TABLE]

Here $e_{i}^{\top}v=v_{i}\in\mathbb{R}$ is a variable that is explicitly kept track of with $v_{i}\approx\mathbb{E}_{k}[g_{k}(W)]_{i}=\sum_{k\neq y_{i}}e^{x_{i}^{\top}(w_{k}-w_{y_{i}}})$ (with exact equality in the limit as $t\to\infty$). Clearly $v_{i}$ in stochastic composition optimization has a similar role as $u_{i}$ has in our formulation for $f$ in (5).

If $i,k$ are sampled with $k\neq y_{i}$ in stochastic composition optimization then the updates are of the form (Wang et al., 2016)

[TABLE]

where $z_{k}$ is a smoothed value of $w_{k}$. These updates have the same numerical instability issues as vanilla SGD on $f$ in (5): it is possible that $\frac{e^{x_{i}^{\top}z_{k}}}{1+v_{i}}\gg 1$ where ideally we should have $0\leq\frac{e^{x_{i}^{\top}z_{k}}}{1+v_{i}}\leq 1$.

Appendix D Proof of general Implicit SGD gradient bound

Proof of Proposition 2.

Let $f(\theta,\xi)$ be $m$-strongly convex for all $\xi$. The vanilla SGD step size is $\eta_{t}\|\nabla f(\theta^{(t)},\xi_{t})\|_{2}$ where $\eta_{t}$ is the learning rate for the $t^{th}$ iteration. The Implicit SGD step size is $\eta_{t}\|\nabla f(\theta^{(t+1)},\xi_{t})\|_{2}$ where $\theta^{(t+1)}$ satisfies $\theta^{(t+1)}=\theta^{(t)}-\eta_{t}\nabla f(\theta^{(t+1)},\xi_{t})$. Rearranging, $\nabla f(\theta^{(t+1)},\xi_{t})=(\theta^{(t)}-\theta^{(t+1)})/\eta_{t}$ and so it must be the case that $\nabla f(\theta^{(t+1)},\xi_{t})^{\top}(\theta^{(t)}-\theta^{(t+1)})=\|\nabla f(\theta^{(t+1)},\xi_{t})\|_{2}\|\theta^{(t)}-\theta^{(t+1)}\|_{2}$.

Our desired result follows:

[TABLE]

where the first inequality is by Cauchy-Schwarz and the second inequality by strong convexity.

∎

Appendix E Update equations for Implicit SGD

In this section we will derive the updates for Implicit SGD. We will first consider the simplest case where only one datapoint $(x_{i},y_{i})$ and a single class is sampled in each iteration. Later we will derive the updates for when multiple classes are sampled, and finally when both multiple classes and multiple datapoints are sampled.

E.1 Single datapoint, single class

Equation (6) for the stochastic gradient with a single datapoint, single sampled class is

[TABLE]

The Implicit SGD update corresponds to finding the variables optimizing

[TABLE]

where $\eta$ is the learning rate and the tilde refers to the value of the old iterate (Toulis et al., 2016, Eq. 6). Since $f_{ik}$ is only a function of $u_{i},w_{k},w_{y_{i}}$, we have the optimal $w_{j}=\tilde{w}_{j}$ for $j\notin\{k,y_{i}\}$ and $u_{j}=\tilde{u}_{j}$ for $j\neq i$. The optimization reduces to

[TABLE]

**Solving for $w_{k},w_{y_{i}}$ with auxiliary variable $b$

**Much of the difficulty in optimizing this equation comes from the interaction between the $e^{x_{i}^{\top}(w_{k}-w_{y_{i}})-u_{i}}$ term and the $\|\cdot\|_{2}^{2}$ terms. To isolate this interaction we introduce an auxiliary variable $b=x_{i}^{\top}(w_{k}-w_{y_{i}})$ and rewrite the optimization problem as

[TABLE]

The inner optimization problem over $w_{k},w_{y_{i}}$ is a quadratic program with linear constraints. Taking the dual and solving yields

[TABLE]

where

[TABLE]

Substituting in the solution for $w_{k},w_{y_{i}}$ and dropping constant terms, the optimization problem reduces to

[TABLE]

We’ll approach this optimization problem by first solving for $b$ as a function of $u_{i}$ and then optimize over $u_{i}$. Once the optimal value of $u_{i}$ has been found, we can calculate the corresponding optimal value of $b$. Finally, substituting $b$ into (10) will give us our updated value of $W$.

**Solving for $b$

**We solve for $b$ by setting its derivative equal to zero in (12)

[TABLE]

Letting $a=x_{i}^{\top}(\frac{\tilde{w}_{k}}{1+\eta\mu\beta_{k}}-\frac{\tilde{w}_{y_{i}}}{1+\eta\mu\beta_{y_{i}}})-b$ and using simple algebra yields

[TABLE]

The solution for $a$ can be written in terms of the principle branch of the Lambert-W function, $P$,

[TABLE]

The optimal value of $b$ given $u_{i}$ is therefore $b(u_{i})=x_{i}^{\top}(\frac{\tilde{w}_{k}}{1+\eta\mu\beta_{k}}-\frac{\tilde{w}_{y_{i}}}{1+\eta\mu\beta_{y_{i}}})-a(u_{i})$.

**Bisection method for $u_{i}$

**Substituting $b(u_{i})=x_{i}^{\top}(\frac{\tilde{w}_{k}}{1+\eta\mu\beta_{k}}-\frac{\tilde{w}_{y_{i}}}{1+\eta\mu\beta_{y_{i}}})-a(u_{i})$ into (12), we now only need minimize over $u_{i}$:

[TABLE]

where and we used the fact that $e^{-P(z)}=P(z)/z$ to simplify the $(K-1)e^{x_{i}^{\top}(\frac{\tilde{w}_{k}}{1+\eta\mu\beta_{k}}-\frac{\tilde{w}_{y_{i}}}{1+\eta\mu\beta_{y_{i}}})-a(u_{i})-u_{i}}$ term. The derivative in (16) with respect to $u_{i}$ is

[TABLE]

where we used the fact that $\partial_{z}P(z)=\frac{P(z)}{z(1+P(z))}$ to work out that $\partial_{u_{i}}a(u_{i})=-\frac{a(u_{i})}{1+a(u_{i})}$.

We can solve for $u_{i}$ using a bisection method. Below we show how to calculate the initial lower and upper bounds of the bisection interval and prove that the size of the interval is bounded (which ensures fast convergence). The initial lower and upper bounds we use depends on the derivative in (17) at $u_{i}=\tilde{u}_{i}$. In deriving the bounds we will use $u_{i}^{\prime}$ to denote the optimal value of $u_{i}$ and $a^{\prime}$ to denote the optimal value of $a$.

**Case: $u_{i}^{\prime}>\tilde{u}_{i}$

**If the derivative is negative then $u_{i}^{\prime}$ is lower bounded by $\tilde{u}_{i}$. An upper bound on $u_{i}^{\prime}$ can be derived from (15):

[TABLE]

where in the inequality we set $a^{\prime}=0$, since the minimal value of $u_{i}$ is monotonically decreasing in $a^{\prime}$. This bound should be used in the bisection method, but for ease of analysis we can weaken the bound:

[TABLE]

where we used the assumption that $u_{i}^{\prime}$ is lower bounded by $\tilde{u}_{i}$ to remove the $(u_{i}-\tilde{u}_{i})^{2}$ term. Thus $\tilde{u}_{i}\leq u_{i}^{\prime}\leq\log(1+(K-1)e^{x_{i}^{\top}(\frac{\tilde{w}_{k}}{1+\eta\mu\beta_{k}}-\frac{\tilde{w}_{y_{i}}}{1+\eta\mu\beta_{y_{i}}})})$. If $(K-1)e^{x_{i}^{\top}(\frac{\tilde{w}_{k}}{1+\eta\mu\beta_{k}}-\frac{\tilde{w}_{y_{i}}}{1+\eta\mu\beta_{y_{i}}})}\leq 1$ then the size of the bounding interval must be less than $\log(2)$, since $\tilde{u}_{i}\geq 0$. Otherwise the gap must be at most $\log(2(K-1)e^{x_{i}^{\top}(\frac{\tilde{w}_{k}}{1+\eta\mu\beta_{k}}-\frac{\tilde{w}_{y_{i}}}{1+\eta\mu\beta_{y_{i}}})})-\tilde{u}_{i}=\log(2(K-1))+x_{i}^{\top}(\frac{\tilde{w}_{k}}{1+\eta\mu\beta_{k}}-\frac{\tilde{w}_{y_{i}}}{1+\eta\mu\beta_{y_{i}}})-\tilde{u}_{i}$. Either way, the size of the interval is upper bounded by $\log(2K)+|x_{i}^{\top}(\frac{\tilde{w}_{k}}{1+\eta\mu\beta_{k}}-\frac{\tilde{w}_{y_{i}}}{1+\eta\mu\beta_{y_{i}}})-\tilde{u}_{i}|$.

**Case: $u_{i}^{\prime}<\tilde{u}_{i}$

**Now let us consider if the derivative in (17) is positive at $u_{i}=\tilde{u}_{i}$. Then $u_{i}^{\prime}$ is upper bounded by $\tilde{u}_{i}$. We can lower bound $u_{i}^{\prime}$ by:

[TABLE]

where the first inequality comes dropping the $(u_{i}-\tilde{u}_{i})^{2}$ term due to the assumption that $u_{i}^{\prime}<\tilde{u}_{i}$ and the second inequality is from the monotonicity of the log function. Recall (13),

[TABLE]

We can upper bound $a^{\prime}$ using the lower bound on $u^{\prime}$:

[TABLE]

Substituting this upper bound for $a^{\prime}$ into (18) and solving yields a lower bound on $u_{i}^{\prime}$,

[TABLE]

Again this bound should be used in the bisection method, but for ease of analysis we can weaken the bound by instead substituting the bound for $a^{\prime}$ into (19) which yields:

[TABLE]

Thus $\log(K-1)+x_{i}^{\top}(\frac{\tilde{w}_{k}}{1+\eta\mu\beta_{k}}-\frac{\tilde{w}_{y_{i}}}{1+\eta\mu\beta_{y_{i}}})-\eta N\gamma_{i}^{-1}\leq u_{i}^{\prime}\leq\tilde{u}_{i}$. The size of the bisection method interval is upper bounded by $\tilde{u}_{i}-x_{i}^{\top}(\frac{\tilde{w}_{k}}{1+\eta\mu\beta_{k}}-\frac{\tilde{w}_{y_{i}}}{1+\eta\mu\beta_{y_{i}}})+\eta N\gamma_{i}^{-1}-\log(K-1)$.

In summary, for both signs of the derivative in (17) at $u_{i}=\tilde{u}_{i}$ we are able to lower and upper bound the optimal value of $u_{i}$ such that interval between the bounds is at most $|\tilde{u}_{i}-x_{i}^{\top}(\frac{\tilde{w}_{k}}{1+\eta\mu\beta_{k}}-\frac{\tilde{w}_{y_{i}}}{1+\eta\mu\beta_{y_{i}}})|+\eta N\gamma_{i}^{-1}+\log(2K)$. This allows us to perform the bisection method where for $\epsilon>0$ level accuracy we require only $\log_{2}(\epsilon^{-1})+\log_{2}(|\tilde{u}_{i}-x_{i}^{\top}(\frac{\tilde{w}_{k}}{1+\eta\mu\beta_{k}}-\frac{\tilde{w}_{y_{i}}}{1+\eta\mu\beta_{y_{i}}})|+\eta N\gamma_{i}^{-1}+\log(2K))$ function evaluations. In practice we use Brent’s method as the optimization routine, which is faster than the simple bisection method. The pseudocode of the entire method is displayed in Algorithm 1.

E.2 Bound on step size

Here we will prove that the step size magnitude of Implicit SGD with a single datapoint and sampled class with respect to $w$ is bounded as ${O(x_{i}^{\top}(\frac{\tilde{w}_{k}}{1+\eta\mu\beta_{k}}-\frac{\tilde{w}_{y_{i}}}{1+\eta\mu\beta_{y_{i}}})-\tilde{u}_{i})}$. We will do so by considering the two cases $u_{i}^{\prime}>\tilde{u}_{i}$ and $u_{i}^{\prime}<\tilde{u}_{i}$ separately, where $u_{i}^{\prime}$ denotes the optimal value of $u_{i}$ in the Implicit SGD update and $\tilde{u}_{i}$ is its value at the previous iterate.

**Case: $u_{i}^{\prime}>\tilde{u}_{i}$

**Let $a^{\prime}$ denote the optimal value of $a$ in the Implicit SGD update. From (14)

[TABLE]

where $u_{i}^{\prime}$ is replace by $\tilde{u}_{i}$ and we have used the monotonicity of the Lambert-W function $P$. Now using the fact that $P(z)=O(\log(z))$,

[TABLE]

**Case: $u_{i}^{\prime}<\tilde{u}_{i}$

**If $u_{i}^{\prime}<\tilde{u}_{i}$ then we can lower bound $a^{\prime}$ from (20) as $a^{\prime}\leq\eta N\gamma_{i}^{-1}$.

**Combining cases

**Putting together the two cases,

[TABLE]

From (10) we know that the step size magnitude is proportional to $a^{\prime}$. Thus the step size magnitude is also $O(x_{i}^{\top}(\frac{\tilde{w}_{k}}{1+\eta\mu\beta_{k}}-\frac{\tilde{w}_{y_{i}}}{1+\eta\mu\beta_{y_{i}}})-\tilde{u}_{i})$.

E.3 Single datapoint, multiple classes.

Consider the case where only one datapoint $i$, but multiple classes $\{k_{j}:k_{y}\neq y_{i}\}_{j=1}^{m}$ are sampled each iteration. Like in Appendix E.1, we will be able to reduce the implicit update to a one-dimensional strongly-convex optimization problem. The resulting problem may be solved using any standard convex optimization method, such as Newton’s method. We do not derive upper and lower bounds for a bisection method as we did in Appendix E.1.

Let us rewrite the double-sum formulation from (5) as $f(u,W)=\mathbb{E}_{i,C_{i}}[f_{i,C_{i}}(u,W)]$ where $i$ is a uniformly sampled datapoint, $C_{i}$ is a set of $m$ uniformly sampled classes from $\{1,...,K\}-\{y_{i}\}$ (without replacement) and

[TABLE]

where $\alpha^{-1}=P(k\in C_{i}|k\neq y_{i})=1-\prod_{j=1}^{m}(1-\frac{1}{K-j})$,

[TABLE]

and $n_{k}=|\{i:y_{i}=k,i=1,...,N\}|$. Using the same derivation as in Appendix E.1, the implicit SGD update is

[TABLE]

The goal is to simplify this multivariate minimization problem into a one-dimensional strongly convex minimization problem. The first trick we will use is to reparameterize $u_{i}=v_{i}-x_{i}^{\top}w_{y_{i}}$ for some $v_{i}\in\mathbb{R}$. This changes the $e^{x_{i}^{\top}(w_{k}-w_{y_{i}})-u_{i}}$ factors to $e^{x_{i}^{\top}w_{k}-v_{i}}$, decoupling $w_{k}$ and $w_{y_{i}}$, which will make the optimization easier. The problem becomes:

[TABLE]

Since $v_{i}=u_{i}+x_{i}^{\top}w_{y_{i}}$ is a linear transformation, (22) is jointly strongly convex in $v_{i}$ and $\{w_{k}\}_{k\in C_{i}\cup\{y_{i}\}}$. Bringing the $w_{k}$ minimizations inside yields

[TABLE]

In Appendix E.1 we were able to reduce the dimensionality of the problem by introducing an auxiliary variable $b$ to separate the exponential terms from the norm terms. We will do a similar thing here. Let us first focus on the inner minimization for $k\in C_{i}$.

[TABLE]

where we have taken the Lagrangian in the final line. The solution for $w_{k}$ in terms of $\lambda$ is

[TABLE]

Thus we know that our optimal $w_{k}$ must satisfy $w_{k}=\frac{\tilde{w}_{k}}{1+\eta\mu\beta_{k}}-a_{k}\frac{x_{i}}{\|x_{i}\|_{2}^{2}}$ for some $a_{k}\in\mathbb{R}$. It can similarly be shown that $w_{y_{i}}=\frac{\tilde{w}_{y_{i}}}{1+\eta\mu\beta_{y_{i}}}+a_{y_{i}}\frac{x_{i}}{\|x_{i}\|_{2}^{2}}$ for some $a_{y_{i}}\in\mathbb{R}$. Substituting this into (23) and dropping constant terms yields

[TABLE]

Using the same techniques as in Appendix E.1 we can analytically solve for the $a$ values:

[TABLE]

where $\sigma(v_{i})=\frac{\eta N\|x_{i}\|_{2}^{2}}{1+\eta\mu\beta_{y_{i}}+\|x_{i}\|_{2}^{2}}\exp\left(\frac{x_{i}^{\top}\tilde{w}_{y_{i}}-v_{i}(1+\eta\mu\beta_{y_{i}})+(\eta N-\tilde{u}_{i})\|x_{i}\|_{2}^{2}}{1+\eta\mu\beta_{y_{i}}+\|x_{i}\|_{2}^{2}}\right)$. Substituting these values into (24) yields

[TABLE]

This is a one-dimensional strongly convex minimization problem in $v_{i}$. The optimal $v_{i}$ can be solved for using any standard convex optimization method, such as Newton’s method. Each iteration in such a method will take $O(m)$ since it is necessary to calculate $a_{k}(v_{i})$, $\partial_{v_{i}}a_{k}(v_{i})$ and $\partial_{v_{i}}^{2}a_{k}(v_{i})$ for all $k\in C_{i}\cup\{y_{i}\}$. The first derivatives are easily calculated,

[TABLE]

as are the second derivatives,

[TABLE]

E.4 Multiple datapoints, multiple classes

Consider the case where $n$ datapoints and $m$ classes are sampled each iteration. Using similar methods to Appendix E.3, we will reduce the implicit update to an $n$ dimensional strongly convex optimization problem.

Let us rewrite the double-sum formulation from (5) as $f(u,W)=\mathbb{E}_{I,C}[f_{I,C}(u,W)]$ where $I$ is a set of $n$ datapoints uniformly sampled from $1,...,N$ (without replacement), $C$ is a set of $m$ uniformly sampled classes from $1,...,K$ (without replacement). The sampled function is of the form

[TABLE]

where

[TABLE]

It will be useful to group the classes that appear in $\cup_{i\in I}\{y_{i}\}$ and those that only appear in $C$:

[TABLE]

The implicit SGD update is

[TABLE]

Like in Appendix E.3, the first step to simplifying this equation is to reparameterize $u_{i}=v_{i}-x_{i}^{\top}w_{y_{i}}$ for some $v_{i}\in\mathbb{R}$ and to bring the $w_{k}$ minimizations inside:

[TABLE]

As done in Appendix E.3, the inner minimizations can be solved analytically by introducing constrained auxiliary variables $b_{ki}=x_{i}^{\top}w_{k}$ and optimizing the dual. We’ll do this separately for $k\in\cup_{i\in I}\{y_{i}\}$ and $k\in C-\cup_{i\in I}\{y_{i}\}$.

For datapoint labels $k\in\cup_{i\in I}\{y_{i}\}$,

[TABLE]

Focusing on the minimization over $w_{k}$:

[TABLE]

The solution for $w_{k}$ in terms of $\lambda_{ki}$ is

[TABLE]

Dropping constant terms, the dual becomes

[TABLE]

where $Q_{k,ij}=\frac{x_{i}^{\top}x_{j}}{1+\eta\mu\beta_{k}}$ and $X_{I}=(x_{i})_{i\in I}\in\mathbb{R}^{D\times n}$ and the optimal $\lambda=Q_{k}^{-1}\left(b_{k}-\frac{X_{I}^{\top}\tilde{w}_{k}}{1+\eta\mu\beta_{k}}\right)$. Now we can solve for $b_{k}$,

[TABLE]

Setting to zero the derivative with respect to $b_{k}\in\mathbb{R}^{n}$ and dividing by 2:

[TABLE]

where $\circ$ denotes the element-wise product, $diag(a)$ is a diagonal matrix, $\mathbf{1}$ denotes the vectors of all ones, $v_{I}=(v_{i})_{i\in I}\in\mathbb{R}^{n}$, likewise for $\tilde{u}_{I}$ and $y_{I}$, and

[TABLE]

Multiplying (27) on the left by $A_{k}^{-1}$, letting $z_{k}=A_{k}^{-1}h_{k}-b_{k}$ and multiplying on the right by $diag(e^{z_{k}})$ yields

[TABLE]

The solution for $z_{k}$ decomposes into separate Lambert-W functions:

[TABLE]

where $P$ is the principle branch of the Lambert-W function applied component-wise. The solution for $b_{k}$ is thus

[TABLE]

where $b_{k}$ is a function of the variable $v_{I}$, which is the only unknown variable that we are yet to minimize over.

For pure class labels $k\in C-\cup_{i\in I}\{y_{i}\}$ the procedure is nearly identical for the the datapoint labels. The optimal value of $w_{k}$ is

[TABLE]

where $\lambda=Q_{k}^{-1}\left(b_{k}-\frac{X_{I}^{\top}\tilde{w}_{k}}{1+\eta\mu\beta_{k}}\right)$ and

[TABLE]

**Final optimization problem

**Substituting the optimal values of $b_{k}$ in (26) yields the final optimization problem

[TABLE]

where $b_{k}(v_{I})$ is from (28). This is a strongly convex optimization problem in $v_{I}\in\mathbb{R}^{n}$. Using standard first order gradient methods, it can be solved to $\epsilon>0$ accuracy in $O(\log(\epsilon^{-1}))$ iterations. The cost per iteration is $O(n^{2}(n+m))$ for the matrix multiplications and $O((n+m)n^{3})$ for the matrix inversions. Note that the matrix inversions do not depend on $v_{I}$ and so they only have to be performed once. Furthermore, if the same minibatches are used each epoch, then the inverted matrices can be calculated just once and stored. The amortized matrix inversion cost is therefore expected to be dominated by the $O(n^{2}(n+m)\log(\epsilon^{-1}))$ cost for solving for $v_{I}$ and the $O(nmD)$ cost of taking the $x_{i}^{\top}\tilde{w}_{k}$ inner products each iteration.

Note that we have assumed that $Q_{k}$ is invertible. As long as the vectors $\{x_{i}\}_{i\in I}$ are independent, this will be the case. If not, then a similar method as above can be developed where a basis of $\{x_{i}\}_{i\in I}$ is used.

Appendix F U-max pseudocode

Appendix G Proof of convergence of U-max method

In this section we will prove the claim made in Proposition 5, that U-max converges to the softmax optimum. Before proving the proposition, we will need a lemma.

Lemma 5**.**

For any $\delta>0$, if $u_{i}\leq\log(1+e^{x_{i}^{\top}(w_{k}-w_{y_{i}})})-\delta$ then setting $u_{i}=\log(1+e^{x_{i}^{\top}(w_{k}-w_{y_{i}})})$ decreases $f(u,W)$ by at least $\delta^{2}/2$.

Proof.

As in Lemma 3, let $\theta=(u^{\top},w_{1}^{\top},...,w_{k}^{\top})\in\mathbb{R}^{N+KD}$. Then setting $u_{i}=\log(1+e^{x_{i}^{\top}(w_{k}-w_{y_{i}})})$ is equivalent to setting $\theta=\theta+\Delta e_{i}$ where $e_{i}$ is the $i^{th}$ canonical basis vector and $\Delta=\log(1+e^{x_{i}^{\top}(w_{k}-w_{y_{i}})})-u_{i}\geq\delta$. By a second order Taylor series expansion

[TABLE]

for some $\lambda\in[0,1]$. Since the optimal value of $u_{i}$ for a given value of $W$ is $u_{i}^{\ast}(W)=\log(1+\sum_{k\neq y_{i}}e^{x_{i}^{\top}(w_{k}-w_{y_{i}})})\geq\log(1+e^{x_{i}^{\top}(w_{k}-w_{y_{i}})})$, we must have $\nabla f(\theta+\Delta e_{i})^{\top}e_{i}\leq 0$. From Lemma 3 we also know that

[TABLE]

Putting in bounds for the gradient and Hessian terms in (29),

[TABLE]

∎

Now we are in a position to prove Proposition 5.

Proof of Proposition 5.

Let $\theta^{(t)}=(u^{(t)},W^{(t)})\in\Theta$ denote the value of the $t^{th}$ iterate. Here $\Theta=\{\theta:\,\|W\|_{2}^{2}\leq B_{W}^{2},u_{i}\leq B_{u}\}$ is a convex set containing the optimal value of $f(\theta)$.

Let $\pi_{i}^{(\delta)}(\theta)$ denote the operation of setting $u_{i}=\log(1+e^{x_{i}^{\top}(w_{k}-w_{y_{i}})})$ if $u_{i}\leq\log(1+e^{x_{i}^{\top}(w_{k}-w_{y_{i}})})-\delta$. If indices $i,k$ are sampled for the stochastic gradient and $u_{i}\leq\log(1+e^{x_{i}^{\top}(w_{k}-w_{y_{i}})})-\delta$, then the value of $f$ at the $t+1^{st}$ iterate is bounded as

[TABLE]

since $\eta_{t}\leq\delta^{2}/(4B_{f}^{2})$ by assumption. Alternatively if $u_{i}\geq\log(1+e^{x_{i}^{\top}(w_{k}-w_{y_{i}})})-\delta$ then

[TABLE]

Either way $f(\theta^{(t+1)})\leq f(\theta^{(t)}-\eta_{t}\nabla f_{ik}(\theta^{(t)}))$. Taking expectations with respect to $i,k$,

[TABLE]

Finally let $P$ denote the projection of $\theta$ onto $\Theta$. Since $\Theta$ is a convex set containing the optimum we have $f(P(\theta))\leq f(\theta)$ for any $\theta$, and so

[TABLE]

which shows that the rate of convergence in expectation of U-max is at least as fast as that of standard SGD.

The proof trivially generalizes to sampling multiple datapoints and classes per iteration by replacing $\log(1+e^{x_{i}^{\top}(w_{k}-w_{y_{i}})})$ with $\log(1+\sum_{j=1}^{m}e^{x_{i}^{\top}(w_{k_{j}}-w_{y_{i}})})$.

∎

Appendix H Results over runtime