Recent neural network and language models have begun to rely on softmax distributions with an extremely large number of categories. In this context calculating the softmax normalizing constant is prohibitively expensive. This has spurred a growing literature of efficiently computable but biased estimates of the softmax. In this paper we present the first two unbiased algorithms for maximizing the softmax likelihood whose work per iteration is independent of the number of classes and datapoints (and does not require extra work at the end of each epoch). We compare our unbiased methods' empirical performance to the state-of-the-art on seven real world datasets, where they comprehensively outperform all competitors.
The paper presents interesting algorithms for minimizing softmax with many classes. The objective function is a multi-class classification problem (using softmax loss) and with linear model. The main idea is to rewrite the obj as double-sum using the dual formulation and then apply SGD to solve it. At each iteration, SGD samples a subset of training samples and labels. The main contribution of this paper is: 1) proposing a U-max trick to improve the numerical stability and 2) proposing an implicit SGD approach. It seems the implicit SGD approach is better in the experimental comparisons. 
I found the paper quite interesting, but meanwhile I have the following comments and questions: 
- As pointed out by the authors, the idea of this formulation and doubly SGD is not new. (Raman et al, 2016) has used a similar trick to derive the double-sum formulation and solved it by doubly SGD. The authors claim that  the algorithm in (Raman et al) has an O(NKD) cost for updating u at the end of each epoch. However, since each epoch requires at least O(NKD) time anyway (sometimes larger, as in Proposition 2), is another O(NKD) a significant bottleneck? Also, since the formulation is similar to (Raman et al., 2016), a comparison is needed. 
- I'm confused by Proposition 1 and 2. In appendix E.1, the formulation of the update is derived, but why we need Newton to get log(1/epsilon) time complexity? I think most first order methods instead of Newton will have linear converge (log(1/epsilon) time)? Also, I guess we are assuming the obj is strongly convex?
- The step size is selected in one dataset and used for all others. This might lead to divergence of other algorithms, since usually step size depends on data. As we can see, OVE, NCE and IS diverges on Wiki-small, which may be fixed if the step size is chosen for each data (in practice we can choose using subsamples for each data). 
- All the comparisons are based on "epochs", but the competing algorithms are quite different and can have very different running time for each epoch. For example, implicit SGD has another iterative solver for each update. Therefore, the timing comparison is needed in this paper to justify that implicit SGD is faster. 
- The claim that "implicit SGD never overshoots the optimum" needs more supports. Is it proved in some previous papers? 
- The presentation can be improved. I think it will be helpful to state the algorithms explicitly in the main paper.