Adaptive Sampling for Efficient Softmax Approximation
Authors: Tavor Baharav
Ryan Kang
Colin Sullivan
Mo Tiwari
Eric Sager Luxenberg
David Tse
Mert Pilanci
Venue: NeurIPS 2024 poster (2024)
Keywords: Multi-armed bandits
adaptive
softmax
attention
Primary area: bandits
TL;DR: Adaptive sampling-based algorithm makes softmax computation sublinear in feature dimension while maintaining PAC guarantees.

Abstract: The softmax function is ubiquitous in machine learning and optimization applications. Computing the full softmax evaluation of a matrix-vector product can be computationally expensive in high-dimensional settings. In many applications, however, it is sufficient to calculate only the top few outputs of the softmax function. In this work, we present an algorithm, dubbed AdaptiveSoftmax, that adaptively computes the top k softmax values more efficiently than the full softmax computation, with probabilistic guarantees. We demonstrate the sample efficiency improvements afforded by AdaptiveSoftmax on real and synthetic data to corroborate our theoretical results. AdaptiveSoftmax yields >10x gain over full softmax computation on most datasets, yielding up to 30x improvement for Mistral7B evaluated on the Wikitext dataset. The adaptive method we propose for estimating the partition function (the softmax denominator) is of independent interest and can be used in other applications such as kernel density estimation.

Adaptive Sampling for Efficient Softmax Approximation
Tavor Z. Baharavt
Ryan Kangt
Eric and Wendy Schmidt Center
Department of Computer Science
Broad Institute
Stanford University Stanford, CA 94305
Cambridge, MA, 02142
baharav@rroaddsstitttee org
txryank@stanfordd edu
Colin Sullivant
Mo Tiwari
Eric Luxenberg Gridmatic Cupertino, CA 95014 eric@gridmatic.o com
AI Division
Department of Computer Science
Software Engineering Institute
Stanford University Stanford, CA, 94305
Pittsburgh, PA 15213
csullivan@sei. cmu.edee
motirari0stanford.oo edu
David Tse Department Electrical Engineering Stanford University Stanford, CA 94305 dntse@stanford.cc edu
Mert Pilanci
Department of Electrical Engineering Stanford University Stanford, CA 94305 pilanci@stanforr edu
Abstract
The softmax function is ubiquitous in machine learning and optimization appli- cations. Computing the full softmax evaluation matrix-vector product can be computationally expensive in high-dimensional settings. In many applications, however, is sufficient to calculate only the top few outputs of the softmax func- tion. In this work, present algorithm, dubbed Adapti SeaptiveSoffi tmax, that adaptively computes the top & softmax values more efficiently than the full softmax computation, with probabilistic guarantees. We demonstrate the sample efficiency improvements afforded by AdaptiveSofte tmax on real and synthetic data to cor- roborate our theoretical results. AdaptiveSof tmax yields > 10x gain softmax computation on most datasets, yielding up to 30x improvement for Mis- tral7B evaluated on the Wikitext dataset. The adaptive method we propose for estimating the partition function (the softmax denominator) is of independent interest and can be used in other applications such kernel density estimation.
full
Introduction
The softmax function appears in many different fields and applications. is often used in multiclass classification problems, as the final operation in a neural network to obtain a probability distribution classes, in reinforcement learning to obtain probability distribution over possible actions, and in statistical mechanics to derive various thermodynamic quantities.
In machine learning applications, the softmax function often appears as the final operation in classification models and in attention layers. Crucially, the softmax function takes a vector weights
+denotes equal contribution
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

as input and returns a probability distribution defined by those weights. Formally, the softmax function for a given temperature parameter defined as:
Su
03 (N) E,eRns
(I)
where μ E IR" the input vector of weights, also referred to as the logits. Usually, the logits are the result of a matrix-vector product (e.g., in fully connected layer where the softmax is used as the nonlinear activation function). The output of the softmax function (Equation 1) is a probability distribution that is "soft" version of the max operator that is differentiable with respect to the logits. The softmax function can thus be used in gradient-based algorithms as a proxy for the non-differentiable max function.
Intuitively, the temperature parameter / controls the peakiness of the softmax output. A larger / corresponds to peakier distribution and "harder" max. The choice of 3 = corresponds to the canonical softmax function 0, and the choice of B oc corresponds to the "hard" argmax. The denominator of Equation (1), Σ; e99,,, is called the partition function and is denoted by Zz.
The softmax function is critical in many popular, recent machine learning applications like large language models (LLMs).. However, it can present computational bottleneck in high-dimensional applications. During the training of neural networks, for example, each training example x requires the computation of the softmax function (x)),, the partition function (z)x), and their gradients. During inference of these models, the number of possible labels for next-token prediction corresponds to the vocabulary size, which can be in the hundreds of thousands for common languages such as English. As such, there has been significant recent interest in accelerating the computation of the softmax function and its derivatives [37, 14, 15].
Key Observations: In many applications, we are only interested in identifying the top few outputs of the softmax function; in these settings, it is unnecessary to compute the smaller entries. This suggests that some of the computation of the full softmax function may be unnecessary and motivates our study, First, we observe that when the input vector the softmax function is the result of matrix-vector product, we can approximate the intermediary computation instead of exactly computing it. This, in turn, allows us approximate the output of the softmax function and converts the problem of computing the softmax function from a computational one to a statistical one. We also note that the softmax output is heavily influenced by the largest input elements which suggests that we can allocate computation adaptively to larger input elements to estimate them with greater certainty. This procedure is inspired by recent work in multi-armes bandits that converts computational problems to statistical ones [4].
Outline: We begin this study with a summary of related work in Section 2. In Section 3, we formalize the reduction of computing the softmax function to a statistical estimation problem. In Section 4, we algorithm based on this reduction. In the same section, we provide
propose the Adapti veeiieeSotttmm probably approximately correct (PAC) guarantees for edaptiveSoftmax and prove that it is more efficient than brute force computation. Crucially, AdaptiveSof resources towards important output values. In Section 5, we demonstrate the empirical advantages of our algorithm in several real-world applications, including in a multiclass classification setting and in large language models. In Section 6, we conclude with a discussion of further applications, potential limitations, and directions for future work.
allocates greater computational
Related Work
Recent work has identified the computational complexity of the softmax function as a significant bottleneck in recent machine learning applications [37]. Well before the attention revolution, [36] proposed methods accelerate softmax computation via hicrarchical model. In their work, a binary tree representing possible outputs as leaves is used and, at each step, the model must predict which path to traverse to a leaf. For a balanced tree with n leaves, the computational complexity of the softmax function is reduced to Ollog n.) from (n)), at the expense of (n)) internal classifiers and providing only an approximate output dictated by the quality of the clustering. Google's 's word2vec models used Huffman trees instead of vanilla binary trees in similar brsed approach [37]. Other approaches include target sampling, noise contrastive estimation and self normalization (summarized [14]). but all of these methods reduce the complexity in terms of the vocabulary size n, rather than

---

by the dimension d. Additionally, our proposed algorithm provides direct PAC guarantees on the original softmax output, instead of approximating the softmax in a sequence of steps without provable accuracy guarantees. Independently, some works have developed fast methods to approximate the softmax gradients during training by using importance sampling over the classes, improving scaling with respect to n [11, 10] leading to a sampled softmax. This is in contrast with Adapti veSof tmax which utilizes importance sampling in an orthogonal direction to subsample the features efficiently, enabling gains in d at both train and test time. These sampled softmax methods were later specialized to kernel-based sampling, resulting in provably bounded bias [12, 38]. However, these and other optimized methods [23] typically require prior knowledge about the desired output label frequencies, leaving them susceptible to phenomena like distribution shift between training and inference data, where the frequency distribution changes at the time the model is evaluated. Unlike these approaches, AdaptiveSe tmax does not require auxiliary knowledge, is adaptive on a per instance basis, and provides provable guarantees for the true softmax computation directly rather than a proxy.
Our algorithm is inspired by adaptive sampling techniques from the multi-armed bandit literature. Randomized algorithms based multi-armed bandit algorithms have seen surge of recent work, due to their ability to provide instance-adaptive guarantees for a variety of problems. This idea was first formalized in the specific case of Monte Carlo Tree Search [26] and later studied in the context of parammer tuning [31]. Recent work has formalized this approach into the framework of Bandit- -Baee Monte Carlo Optimization [4]. where the computational task is reduced to of statistical estimation that is solved efficiently with adaptivity. Applications of this framework include finding the medoid of dataset [5, 7, 41], k-nearest neighbor graph construction [30, 34), Monte Carlo permutation- hased multiple testing [46], and an adaptive singular value decomposition (SVD) [24). Most relevant is the recent work of [8], where the authors provide general framework for adaptive sampling to approximate function evaluation at unknown but estimable points. This work provides general guarantees, but requires bound on the Lipschitz factor of the function's 's gradients as input and has potentially poor performance on specific function classes due to its generality.
A sub-problent in our softmax approximation is identifying the index of the largest component; this is equivalent to the Maximum Inner Product Search (MIPS) problem on the preceding matrix- vector product. MIPS is common problem that has inspired many directions of research [22, 32]. Many of these algorithms focus on specific use cases and place restrictive assumptions on the data (c.g., that all elements of the matrix- vector product underlying data distribution, or lack PAC guarantees. One large family of MIPS algorithms are based hashing (LSH) [19, 2]. In addition to significant preprocessing overhead and practical implementation issues, a shortcoming of these LSH-based approaches is that the maximum dot product is often small compared to the vector norms in high dimensions, which necessitates many hashes and significant storage space (often orders of magnitude more than the data itself). Promising LSH-based algorithms have recently been applied to the problem of softmax computation [1, 15]. These methods focus on intensive preprocessing and work primarily by attaining gains in terms of n. In contrast, AdaptiveSof tmax subsamples matrix- -ector products and obtains gains with respect to d. Furthermore, Adapti ve provides an instance-adaptive algorithm with no required preprocessing and still has PAC guarantees.
positive), require preprocessing. are not adaptive the
on locality-sensitiee
3 Problem Formulation
In this work, we focus on the problem of identifying the k: largest entries of the softmax of an input vector that is the result of a computationally expensive matrix-vector multiplication. Specifically, we analyze the setting where the input vector j is the result of matrix vector product Axr, as is common in the final linear layer of neural networks (such scenarios frequently arise in other machine learning problems well [4]). Our objective is to design an algorithm that can with probability at least 1 - d estimate the top k values to multiplicative accuracy e, where e and 8 are given input parameters. For clarity of exposition, we focus on the case of k = 1, i.e., identifying and estimating the largest component. All our theoretical results, however, easily extend to the setting & > (discussed in Section 4.2).
Notation:
We use [n] to denote the set otherwise specified. We use ||, to denote the Orlicz norm (i.e., the vub-Gaussianity parameter or variance proxy) of a random variable; this is discussed in greater detail in Appendix A.2 [44]. For matrix A and vector z, we denote the resulting product as μ Az. Assuming for notational
nn and ||| to denote the vector &2 norm, unless
3

---

simplicity that the arms are in sorted order > A2 2 ... ... we define the gaps between the entries of μ as A; = μi p4. We use the convention from the best-arm identification literature that A 22 and assume that A2 > O (this assumption is easily relaxed). Furthermore, we define a; = e and % = @ which are proportional to the optimal first (respectively, second) order sampling frequencies; these are discussed further in Section B.2. Finally, we define p as the softmax output, and i" as its largest entry (assumed to be unique), i.e.,
(A(Ax) P, i" argmaxp:.
*E|n]
Our goal is design an algorithm which efficiently outputs the best index i" and an estimate its value where, with probability at least 1 - 6, the best index is correct and the estimated value is within factor of multiplicative accuracy. Mathematically. defining the algorithm's 's outputs as d [n]] and P.- E [0, 1], we define the success events Ex, Eest where the algorithm identifies the largest entry, and where it estimates its value within multiplicative accuracy We define the algorithm as providing (s, 6)-PAC guarantees if these events happen simultaneously with probability at least 1 6, with respect to the randomness of the algorithm.
Eu = f =
Fest = {(1 E)p ft- S ( e)p:- IP (Bia n Ees) 2 5
(3)
Our objective them becomes to design algorithm that satisfies Equation (4) and minimizes the requisite sample complexity.
4 AdaptiveSoftme Algorithm
We now introduce the idaptiveSooftaxx.llaeae max in Algorithm First, AdaptiveSof approximates the softmax normalization constant Za to a multiplicative accuracy of e/4 via AdaptiveSof
which approximates the output of the soft-
(Algorithm 2). Next, identifies the best arm (or top & arms, depending on the setting) using standard multi-armed bandit algorithm, BestArmId (Algorithm 3). In our setting, "arms" correspond to computing coordinate- wise scalar product
different rows of and pulling arm corresponds Ai. for some coordinate (we provide more formal overview of the best-arm identification problem and the associated algorithm in Appendix A). Finally, Adapti iveof tmax estimates the value of the identified best arm (or top & arms) to a multiplicative accuracy of e/4 by sampling each arm a sufficient number of times via Est EstimateArm (Algorithm 10).
We prove (s, 5) -PAC guarantees for AdaptiveSoftma by mmion-bounding over the error prob- abilities in each step of Algorithm Our results will show that, with probability least AdaptiveSor tmax is able to identify the largest output of the softmax function and estimate its value multiplicative accuracy c.
Algorithm 1 Adaptive Sof tmax
Input: Matrix A, vector z, temperature 3, error c, failure probability variance proxy o²
2: Output: P.. and f², highest softmax probability and its index
3: # Estimate denominator of softmax
4: Ż
5: # Compute index of best arm.
6: & --BestArmla (A, (A, A, 5/3,
7: # Estimate value best arm
8: f.z EEttimateArm A,,, x,, e,,,, 5/355/3)
9: p:- = cxp(
10: return 6-
These inputs are typical in the multi-armed bandit setting, but the variance proxy x² merits additional discussion. In order for our random- random-sampling samp-ing-mmllaeee approach to succeed, a bound the rate of

---

Algorithm 2 Normaal
1: Input: Matrix A, vector æ, temperature 3, target error e, failure probability variance proxy o²
2: Output: Zg. estimate the partition function
3: Compute üú; using To 1732d² (6n/66) coordinate samples for each arm
2x lood sD
4: C f-{)
5: 0; f- (
(()-
6: %/ f- --
7: Sample each arm n, = min(n;, d) times to recompute the estimates fi. where
6m 16 log
9s) 2%
π4 = 022 max 17
6n
8
16log 0
12
-2 âi
Σg āj
Ej Ejas)
8: return ZB = E, e
concentration of the estimators μå; is required; the quantity o² governs the concentration rate, discuss in Appendix In practice, such bound holds very generally, for example as long as A and have bounded entries. For algorithmic simplicity utilize the following assumption.
Assumption 1. We assume that we are given variance proxy bound o² for the sub- SabbGaassian parameters of the constructed estimators:
> ||AraaJ Vi, for J ~ Unif((n]).
We provide theoretical guarantees for AdaptiveSoftnas under Assumption 1. Recall that we defined our optimal first and second order sampling frequencies og = ehd. and " @/u/22 (see Appendix B.2). We first show in Proposition that softmax normalization estimation algorithm (Algorithm 2) obtains the desired guarantees.
Proposition For input E E 1/,) E (0, 1). and satisfying Assumption 1, Algorithm will, with probability least estimate 73 E3 eavo, multiplicative accuracy of c. On this
success event, Algorithm 2 requires at most
2
CB8 nog (7)
+ log
Σni 73
Σ
oj
+ log
samples for some absolute constant C, where non-asymptotic bounds with numerical constants are provided in Appendix B.
With the sample complexity of Algorithm 2 bounded, the complexity of best arm identification and the cost of estimating the best arm to target accuracy are readily available from the multi-armed bandit literature. This enables to state overall result for Adapti seetriiieoottmm tmax (Algorithm 1) in the following Theorem.
Theorem For input E (0, 1/2), 8 (0. 1), and a satisfying Assumption I, Algorithm identifies the largest component in 0 (Az) and estimates its value multiplicatiie accuracy of with probability at least 1 8. as in (4). On this success event. the algorithm uses samples where
log
log d
B log (3)
Ba log(1/8)
Co2 C
T T
11 log
+
A?
E, 0g
e2
for some absolute constant C. Tighter bounds with non-asymptotic numerical constants are provided in Appendix B.
The proofs of these two results are detailed in Appendix B; we provide some intuition and brief sketches of the proofs here.
5

---

For Proposition 1, we first show that we can estimate the quantities {az}, {ri) to constant multi- plicative error with high probability. This allows to construct a sampling scheme based off of the asymptotically optimal sampling frequencies, and guarantee that each arm is sampled at least half of what this asymptotically optimal frequency requires. Then, sampling each arm enough so that the first order Taylor expansion of e is sufficiently accurate, we can sample further according to these determined frequencies to guarantee PAC estimation. This is an improved and specialized modification of [8] that exploits the structure of the softmax function to remove the assumption of Lipschitz gradients and yield improved sample complexity (this is discussed further Appendix B.3).
Next, we utilize a classical best-arm identification algorithm to identify the best arm with high probability, leveraging standard results in Bandit-Based Monte Carlo Optimization [6]. Finally, we sample the identified best arm enough times to estimate its value to multiplicative accuracy e/4 with high probability. By union bounding over these error probabilities, we achieve the desired PAC guarantees.
4.1 Interpreting Theoretical Results
We now simplify and further interpret the sample complexity results in Theorem First, note that the E dependence exhibited by
on (Algorithm 2) is optimal: it is inherent even in estimating the mean of the best arm to accuracy E. The cost stemming from the second order error, which scales as e e-1, is bounded between 8820² log(n/8)s" and n.8² log(n/5)e where in the case where one arm is much better than the rest this will match the first term. Concretely, we analyze the setting where the minimum gap is A, i.e. P4 μ4 A, 2 A for all i.
Corollary 1. Under the conditions of Theorem I, when the minimum gap is at least A, Algorithm 1 identifies and provides (e, 8 estimation (Equation (4)) of the largest softmax entry, using
T T C 8 log
T1 8
n 1n1 n + eB2 + =f + no log
n2 lologe d
A
A
samples for some universal constant C. In the case where the gap is large (A 2 3 logn), B is not too small, and e (see Equation (47) for precise statement), this can be simplified to
TSC (oo (0)
n (n + + +e 2lg((181)
where all sample complexities are for the 1 success event.
The proof of this upper bound is in Appendix B.I. Note that this directly implies that when the gap is large (i.e. there is a clear largest output element) and E is constant, the sample complexity is nearly linear in and is upper bounded by C log(n/O)
4.2 Implementation details and extensions
There are many techniques that we can use to extend and improve AdaptiveSoftmax in practice. We discuss changes from the written algorithm in detail in Appendix C.
Randomized Hadamard Transformation: The variance- peroice-poxyyoao bound 0² of the arms plays large factor in the AdaptiveSof tmax algorithm's 's sample complexity. The underlying sub- Sab-Gaussiantty parameter of these estimators can be improved using techniques from randomized numerical linear algebra, such as the randomized Hadamard transform [42]. If a small number of entries dramatically increase the variance of the estimator, then the randomized Hadamard transform will make the coordinates more uniform. We provide theoretical guarantees for this approach in Appendix A.3.1.
Top-k Identification: Extending our algorithmic results from best arm identification (top 1) to identifying multiple components (top k) follows directly from existing multi-armed bandit algorithms. Numerous algorithms have been developed for this setting [21], and variants for computational settings have been developed and studied in [6]. For simplicity and clarity, we focused on the top I identification this paper, but the top à extension readily follows. Furthermore, in numerical experiments we observe estimating the normalization constant Z3 dominates the sample complexity. and the increase in cost from identifying the top arms and estimating their values to multiplicative accuracy e/ 4 is minimal.

---

Relaxing Assumption of Known sub-Gaussian Parameter o²: Assumptions regarding known arm concentration parameters are common multi-armed bandit works and simplify theoretical exposition. These results can naturally be extended in several directions. One simple extension is to the setting where we have a separate Gabbsaaassa parameter of for each arm, i.e., heterogeneous variances. A more practical extension is to the setting where we do not have bound on the sub- Gaussianity parameter for each arm but know that the arm returns are bounded. this setting, common multi-armed bandit approach is to utilize the empirical variance [35]. These approaches are discussed further [8].
Improved Estimators ji: Naively, the Adapti iveof tmax algorithm samples coordinates uniformly random with replacement from the of coordinates {1..... d} estimate each Σ, At This procedure can be improved in several ways. For example, we may utilize importance sampling and sample each coordinate with probability 2j & ææy|. Furthermore, we can sample coordinates without replacement; this known to yield tighter confidence intervals than sampling with replacement [9]. We can combine these techniques and compute the effective variance as in [18]. Sampling without replacement can be achieved in computationally efficient manner via Gumbel sampling [27]. We discuss these details further in Appendix A.3; these details may be of independent interest.
Experiments
In this section, we demonstrate the empirical advantages of AdaptiveSoftman over the brute-force softmax computation in terms of sample complexity. All of our results are reproducible via I-line script, publicly available on GitHub github. comm Twtrrhrrruuuupp///o adaptiveSof torrex...aa...ommm
5.1 Complexity on Synthetic Data
Crucially, the AdaptiveSof tmax algorithm scales sublinearly in d. More precisely, Corollary implies that, for fixed and 6, the sample complexity of the AdaptiveSoftmax algorithm scales as Oog n). We now empirically validate this behavior.
We first run the ddaptiveSoftmat algorithm on two synthetic datasets. In each dataset, we generate A and æ with 100 and vary d.
In the first synthetic dataset, we set to be a d-dimensional vector of all 1s. We draw each element of A iid. (O 1) and add the vector of all the first row of A, thereby planting signal. In expectation, the first row of A will have inner product d with æ whereas all other rows will have inner product 0 with æ. Furthermore, all arms have expected variance of that scales with d.
naive
Sample complexity
AclaptiveSoftmax
Sample complexity
10a
naive AdsptiveSortmax
10a
10a
20k
30k
40k
50k
60k
70k
80k
30k
50k
60k
80k
90k
Dimension
Dimension
(a) Scaling Baseline: All-Ones Query
(b) Scaling Baseline: Sign Query
Figure 1: Sample complexity of the AdaptiveSoftmex algorithm and the brute-force softmax computation on two different synthetic datasets as a function d. Error bars are obtained from 100 random trials. The sample complexity of the AdaptiveSof tmax algorithm scales with respect to d for (a) but does not for (b), as expected. The average gains for 10% and æ = 30% are 3.953x for (a) and 29.039x for (b), increasing with increasing dimension. Confidence intervals are 1std.
In the second synthetic dataset, draw each element of A Lid \/ and set be |A.. the entrywise absolute value of the first row of A. Here, arms have expected variance of O

---

Figures (() and 1(b) demonstrates the scaling of the AdaptiveSoftmax Algorithm on each of the two datasets. On the first synthetic dataset, the AdaptiveSoftmax algorithm scales with d because the variance proxies of do. On the second synthetic dataset, the Ad tmax algorithm does not exhibit significant scaling with d. On both datasets, the AdaptiveSof tmax algorithm significantly outperforms the naive brute-force computation of the softmax function.
5.2 Multinomial Logistic Regression
Multinomial logistic regression (MNL) is form of multiclass classification in which the final operation performed by the classifier is of the form:
B(w ((W))()) h(o))
P(y = c)- Ee=1 ((nT(()))
(5)
i..,, the probabilities that datapoint belongs each class e is given by the softmax applied the vector (h()), where W is the matrix containing rows w1, We and h(x) is a latent representation of x (i.e., the forward pass of some neural network on r))
The multinomial logistic regression naturally amenable to accelerated softmax computation in Equation (5). In many real-world settings, both the number of classes C' and the dimension of the latent representation h(x) (and therefore the dimensionality of each w.) can very large, motivating the usage of AdaptiveSof tmax to identify and estimate the probability of the most likely class. However, the application of Adapti AdaptiveSof tmax extends far past vanilla MNL. For instance, the final layer of any neural network classifier utilizing softmax also be viewed now provide several such practical settings for which we demonstrate the benefits of applying the AdaptiveSof tmax algorithm.
MNL problem. We
5.3 AdaptiveSoftmax Performance on Real Data
We now demonstrate the performance of the AdaptiveSoftmax algorithm on several real-world datasets. For each setting, we provide the sample complexity gain relative to the sample complexity of the brute-force, naive softmax computation sample complexity nd. We also provide the success rate of our algorithm in each setting, i..., the proportion of times the adaptiveSof tmax algorithm correctly identifies the maximum likelihood output (i.e. 8 =i+) and estimates its probability pe- within a multiplicative error of c = 30%.
5.3.1 Application to CNNs
We consider the application of AdaptiveSoftmax to CNN classifiers two distinct image classifi- cation datasets:
1. The MNIST dataset, containing black and white images of handwritten digits as input and ten output classes representing all ten possible digits.
2. The EuroSAT dataset, containing RGB satellite imagery as input and ten output classes, representing possible land types (e.g., river, residential. etc)
On both of these datasets and for distinct architectures, we show that AdaptiveSoftman provides drastic improvement in sample efficiency.
For the MNIST dataset, we train shallow CNN from scratch with two convolutional blocks (Conv2d, ReLu, MaxPool, BatchNorm). This model achieves over 99% accuracy on the test set. The matrix is obtained by extracting the weight matrix of the model's final linear layer. The vector is extracted as the output of the final hidden layer (the layer before the final linear layer) constructed by passing the MNIST image through the trained model and flattening the result. The dimensionality of is adjusted by changing the number of output channels of the convolution blocks. The sample complexity of our algorithm is measured by running the algorithm on 1000 different images in test set with same matrix A. The empirical error rate d is calculated as the fraction of experiments where the adaptive algorithm fails to identify the same class, or fails to estimate the probability to accuracy e, as assigned by exact computation.
MNIST

---

EuroSAT We also utilize a larger pre-trained CNN classifier fine-tuned on the EuroSAT dataset, to show that AdaptiveSof tmax works with larger more sophisticated CNNs. Specifically, we freeze all convolution blocks of VGG-19 (pretrained on ImageNet) and changed the final output dimension to 10 classes for FuroSAT without freezing the weights. The resulting model achieves 92% accuracy on the test set. As before, the matrix A can be extracted from the weights of the final linear layer and the vector z represents the final hidden layer activations. The empirical error rate is calculated in the same manner as for MNIST.
Dataset (Model)
8 = 10%
ò = 5%
0 = 1%
EuroSAT (VGG-19)
5.18x (80.62%)
5. 5.16x (83.00%
4.54x (98. (98.7%%)
MNIST (Shallow CNN)
8.95x (92.25%)
8.81x (93.75%)
8.13x (99.38%)
Table 1: Performance improvement and success rate afforded by Adapti veeof tmax for multinomial logistic regression on two different real-world datasets. We used a total of q = 800 test queries to measure success rate.
5.3.2 Application to LLMs
We also apply the AdaptiveSoftmax algorithm to LLMs using HuggingFace's AutoModee1FoCCusall.. module for the task-generation task [45]. The matrix A the 1m- in-head layer for each model, and the queries z are the final hidden states of the model that is extracted by running a forward pass of the model on the given dataset with a window moving at a certain stride. The context window and stride is modified to generate desired number of queries.
Dataset (Model)
= 10%
= 5%
8 = 1%
Wikitext (GPT-2)
8.25x (88.94%) 9
7.80x (93.54%)
6.67x (98. (6%)
Wikitext (Llama3-7B)
14.68x (91 )
11.43x (94.04%)
6.88x (99.38%
Wikitext (Mistral7B)
32.65x (89.08%)
26.37x (91.20%)
17.71x (97.77%
Penn Treebank (GPT-2)
8. 10x (81.68%)
7.50x (90.73%)
6.66x (96.7 (96.9%%
Penn Treebank (laama377B)
19.18x (87.82%)
16.57x (91.60%)
10.72x (97.81%)
Table 2: Performance improvement and success rate afforded by AdaptiveSoftmax for LLM inference (improvement for final softmax layer). Experiment details in Section 5.3.2. We used q = 1000 unseen test queries to measure d-accuracy.
Our matrix A is the extracted Im-head from HuggingFace's AutoModelForCCusaalL.M for the four models: GPT-2 (n 50257, d 768). Llama3-7B (n 128256, d 4096), Mistral7B (n = 32000, d = 4096), and Gemma7B (n = 256000, d = 3072). Our task is fask-generation, and we generate our queries x by using two datasets (Wikitext and Penn Treebank) with a sliding window of certain stride. Stride and context window is set to get y = 1000 number of queries. Constants and confidence intervals given by theory empirically quite loose, so we tuned algorithm parameters (constant coefficients for stage length and confidence interval width) on initial training data, described in Appendix C. An aggressive tuning strategy was undertaken order demonstrate the potential gains in sample complexity provided by AdaptiveSof tmax Specifically, constant multiples were applied to variance estimate within Algorithm and Algorithm 2. Due the limited sample set, this approach occasionally overoptimized the constants on training data, yielding lower success rates than targeted. However, from the results, it clear that this target parameter still provides users sufficient control over the tradeoff between true success rate and sample complexity.
Discussion, Limitations, and Future Work
In this work, we proposed theoretically novel and practically efficient algorithm for approximating the softmax function. We provided theoretical guarantees on the accuracy of our approximation and demonstrated that, with fewer samples than exact computation, we approximate the softmax function to the desired accuracy. We further demonstrated the viability of our proposed algorithm in two real-world settings, multinomial logistic regression and LLM inference.
potential limitation of our proposed algorithm is that it is most beneficial when the inner dimension of the matrix vector product is high dimensional; its benefits over exact computation are more modest

---

when the inner dimension is small. In particular, the exact computation of the matrix-vector product preceding a softmax operation is usually performed efficiently using BLAS (Basic Linear Algebra Subroutines, which are highly optimized). Adaptivity at its core is inherently sequential, whereas BLAS operations take advantage of batch computation. algorithms, with only a logarithmic number of rounds of adaptivity, but there are important directions of future work to realize these thcoretical gains in practice.
this work we proposed minimally adaptive
Limitations: Theoretical sample complexity bounds are useful for understanding the fundamental properties of an algorithm, but in practice, wall-clock time is often the metric of interest. Many of the steps in our algorithm can be batched and made BLAS efficient, yielding comparable wall clock times brute force computation. However, in general adaptivity is the opposite of batching, as can be seen when we modify our algorithm to adapt to arm specific variances. In this case, we must sample each arm individually, as the number of samples required for each arm is different. This is trade-off between adaptivity and wall-clock time, and in practice, the choice of which to prioritize depends on the specific application (energy efficiency, computational resources, etc.). There are also possible theoretical analyses, where we e.g. create batches of arms with similar empirical variance and sample all arms within a batch together, leading to trade-off between adaptivity and batched computational efficiency. Additionally, in large language models, the final softmax layer is often not a computationally significant step, so while such a method may greatly accelerate multinomial logistic regression, more work may be required to have this accelerate LLMs.
Given the ubiquity of the softmax function in today's 's machine learning workflows, we hope that our algorithm will help pave the way for an optimized adaptive softmax that can accelerate wide class of machine learning models. An interesting direction of future work is trying to combine this multi-armed bandit approach with LSH [15] to obtain (for the attention case) subquadratic complexity in n, and sublinear complexity in d. The adaptive method we propose for estimating the normalizing constant of the softmax function is of independent interest, and holds potential for applications in kernel density estimation and other machine learning tasks.
Acknowledgements
Mert Pilanci was supported in part by the National Science Foundation (NSF) under Grant DMS- 2134248; in part by the NSF CAREER Award under Grant CCF-2236829; in part by the U.S. Army Research Office Early Career Award under Grant 1911NF-21-1-0242 and in part by the Office of Naval Research under Grant 1--111---4-2144, Tavor Baharav was supported by funding from the Eric and Wendy Schmidt Center at the Broad Institute of MIT and Harvard.
References
[1] Josh Alman and Zhao Song. "Fast attention requires bounded entries". In: Advances in Neural Information Processing Systems 36 (2024).
[2] Alexandr Andoni et al. "Practical and optimal LSH for angular distance". In: Advances in neural information processing systems 28 (2015).
[3] Jean-Yves Audibert, Rémi Munos, and Csaba Szepesvári. "Exploration cxxlotttonnnnlloottaiiiinno tradeoff using variance estimates in multi-armed bandits". In: Theoretical Computer Science 410.19 (2009), pp. 1876-1902.
[4] Vivek Bagaria et al. "Bandit-Based Monte Carlo Optimization for Nearest Neighbors". In: IEEE Journal on Selected Areas Information Theory (2021).
[5] Vivek Bagaria et al. "Medoids in almost-linens time via multi-arman bandits". In: International Conference on Artificial Intelligence and Statistics (2018), pp. 500-509.
[6] Tavor Baharav and Tze Leung Lai. "Adaptive Data Depth via Multi-Armed Bandits", In: Journal of Machine Learning Research 24. 24.155 (2023), pp. 1-29.
[7] Tavor Baharay and David Tse. "Ultra fast medoid identification via correlated sequential halving". In: Advances Neural Information Processing Systems 32 (2019).
[8] Tavor Baharav et al. "Approximate Function Evaluation via Multi-Armed Bandits". In: Inter- national Conference on Artificial Intelligence and Statistics. PMLR. 2022, pp. 108-135.
[9] Rémi Bardenet and Odalric- Maillard. "Concentration inequalities for sampling without replacement". In: Bernoulli 21.3 (2015), pp. 1361-1385.
10

---

[10] Yoshua Bengio and Jean- Sén--éébsttinn Senécal. "Adaptive importance sampling to accelerate training of a neural probabilistic language model", In: IEEE Transactions on Neural Networks 19.4 (2008), pp- 713-722.
[11] Yoshua Bengio and lean-Sébastien Senécal. "Quick training of probabilistic neural nets by importance sampling". In: International Workshop on Artificial Intelligence and Statistics.
PMLR. 2003, pp. 17-24.
[12] Guy Blanc and Steffen Rendle. "Adaptive sampled softmax with kernel based sampling". In: International conference on machine learning. PMLR. 2018, pp. 590-599.
[13] Sébastien Bubeck, Nicolo Cesa-Bianchi, et al. "Regret analysis of stochastic and nonstochastic multi- bandit problems". In: Foundations and Trends@ in Machine Learning 5.1 (2012),
Pp. 1-122.
[14] Wenlin Chen, David Grangier, and Michael Auli. "Strategies for Training Large Vocabulary Neural Language Models". In: Proceedings of the 54th Annual Meeting of the Association for Computational Linguistics. 2016, pp. 1975-1985.
[15] Insu Han et al. "HyperAttention: Long-context Attention in Near-Linear Time". In: arXiv preprint 1 (2023).
[16] Eshcar Hillel et al. "Distributed exploration in multi-armeds bandits". In: Advances in Neural Information Processing Systems 26 (2013).
[17] Ari Holtzman et al. "The curious case of neural text degeneration". In: arXiv preprint arxiv: 097.09751 (2019).
[18] Daniel G Horvitz and Donovan J Thompson. "A generalization of sampling without replace- ment from finite universe". In: Journal of the American statistical Association 47.260 (1952),
pp. 663-685.
[19] Piotr Indyk and Rajeev Motwani.
"Approximate nearest neighbors: towards removing the curse of dimensionality". In: Proceedings of the thirtieth annual ACM symposium on Theory of computing. 1998. pp. 604-613.
[20] Andrei Ivanov et al. "Data movement is all you need: A case study on optimizing transformers".
In: Proceedings of Machine Learning and Systems (2021), pp. 711-732.
[21] Kevin Jamieson and Robert Nowak. "Best-arm identification algorithms for multi-armed bandits in the fixed confidence setting". In: Annual Conference on Information Sciences and Systems. 2014, pp. 1-6.
[22] Hervé Jégou et al. "Searching in one billion vectors: re-rank with source coding". In: 2011 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE. 2011, pp. 861-864.
[23] Armand Joulin et al. "Efficient softmax approximation for GPUs". In: International conference on machine learning. PMLR. 2017, pp. 1302-1310.
[24] Govinda Kamath, Tavor Baharav, and Ilan Shomorony. "Adaptive learning of rank-one models for efficient pairwise sequence alignment". In: Advances in Neural Information Processing Systems 33 (2020), pp. 7513-7525.
[25] Zohar Karnin, Tomer Koren, and Oren Somekh. "Almost optimal exploration in multi-armed bandits".. In: International Conference on Machine Leurning. 2013, pp. 1238-1246.
[26] Levente Kocsis and Csaba Szepesvári. "Bandit based monte-carlo planning". In: European
conference on machine learning. Springer. 2006, pp. 282-293.
[27] Wouter Kool, Herke Van Hoof, and Max Welling. "Stochastic beams and where to find them: The gumbel-top-k trick for sampling sequences without replacement". In: International Conference on Machine Learning. PMLR. 2019, pp. -409-35508
[28] Tze Leung Lai and Herbert Robbins "Asymptotically efficient adaptive allocation rules". In: Advances in applied mathematics 6.1 (1985), pp. 4-22.
[29] Tor Lattimore and Csaba Szepesvári. Bandit algorithms. Cambridge University Press, 2020.
[30] Daniel LeJeune, Reinhard Heckel, and Richard Baraniuk. "Adaptive estimation for approxi- mate k-nearest-neighhor computations". In: The 22nd International Conference Artificial Intelligence and Statistics. PMLR. 2019, pp. 3099-3107.
[31] Lisha Li et al. "Hyperband: A novel bandit-based approach to hyperparameter optimization". In: The Journal of Machine Learning Research 18.1 (2017), pp. 6765-6816.
11

---

[32] Stephan S Lorenzen and Ninh Pham. "Revisiting wedge sampling for budgeted maximum inner product search". In: Joint European Conference on Machine Learning and Knowledge Discovery in Databases. Springer. 2020, pp. 439-455.
[33] Pinyan Lu, Chao Tao, and Xiaojin Zhang. "Variance-depenee dependent best arm identification". In: Uncertainty in Artificial Intelligence. PMLR. 2021, pp. 1120-1129.
[34] Blake Mason, Ardhendu Tripathy. and Robert Nowak. "Nearest neighbor search under uncer- tainty". In: Uncertainty in Artificial Intelligence. PMLR. 2021, pp. 1777-1786.
[35] Andreas Maurer and Massimiliano Pontil. "Empirical bernstein bounds and sample variance penalization", In: arXiv preprint arXiv:0907.3740 (2009).
[36] Frederic Morin and Yoshua Bengio. "Hierarchical probabilistic neural network language model".. In: International workshop on artificial intelligence and statistics. PMLR. 2005,
pp. 246-252.
[37] Kezban Dilck Onal et "Neural information retrieval: At the end of the early years". In: Information Retrieval Journal 21 (2018), pp. 111-182.
[38] Ankit Singh Rawat et al. "Sampled softmax with random fourier features". In: Advances in Neural Information Processing Systems 32 (2019).
[39] Max Simchowitz, Kevin Jamieson, and Benjamin Recht. "The simulator: Understanding adaptive sampling in the moderate-confidence regime". In: Conference on Learning Theory.
PMLR. 2017, pp. 1794-1834.
[40] Aleksandrs Slivkins et al. "Introduction to multi-armed bandits". In: Foundations and Trends@ in Machine Learning 12.1-2 (2019), pp. 1-286.
[41] Mo Tiwari et al. "BanditPAM: Almost linear time k-medoids clustering via multi-armed bandits". In: Advances Neural Information Processing Systems 33 (2020), pp. 10211-10224
[42] Joel A Tropp. "Improved analysis of the subsampled randomized Hadamard transform". In: Advances in Adaptive Data Analysis 3 (2011), pp- 115-126.
[43] Tim Vieira. Gumbel-max trick and weighted reservoir sampling. 2014. URL: http / timvieira. github io / blog / post / 2014 / 08 / 01 gumbel max trick and reighttedrreseevirr-ssee swwporrrr.aaann/.........m..
[44] Martin J Wainwright. /iigh-dimensional statistics: A non-asymptotic viewpoint. Vol. 48. Cam- bridge university press, 2019.
[45] Thomas Wolf et al. "Huggingface's transformers: State-of-the-art natural language processing".
In: arXiv preprint arXiv: 1910.0377 (2019).
[46] Martin Zhang, James Zou, and David Tse. "Adaptive monte carlo multiple testing via multi-
armed bandits". In: International Conference Machine Learning. Proceedings of Machine Learning Research. 2019, pp. 7512-7522.
12

---

A
Bandit preliminaries
To make this work accessible to a broad audience, we provide self contained introduction to the multi-armed bandit setting.
A.1 Best arm identification
We consider stochastic multi-armed bandit problem [13, 40, 29] with n arms (distributions). where each arm i has an unknown mean reward μ4, At each time step t. the algorithm selects an arm 4 [n] and receives reward X,, drawn from the distribution of arm 4t Early work in the multi-armed bandit literature focused on the regret minimization setting. where the goal is to maximize the cumulative reward (sum of arm pulls observed so far), motivated by applications such as online advertising and gambling [28]. Recent work has seen increased interest in the best-arm identification setting, where the goal to identify the arm with the highest mean reward with high probability, motivated by applications such clinical trials. Significant research has been devoted to obtaining optimal logarithmic factors, but for the sake of clarity we highlight here a simpler and empirically well performing algorithm, multi-roume -ramm from [25].
Algorithm 3 BestArmId (modification of best-arm identification from [25])
Input: n arms, error probability 6, variance proxy o²
Output: Best arm i* with probability at least -
So + [n]
to f-
while | do
r+r+1 G 2-
f,
+ 1 log(lorg442///
for all arms E S, do
Pull arm 4, - 4r- times and observe rewards Xi,, 14by Xi
Per t 1 E X1,
b Update mean estimates
f,
Cirr / 20 log(4nr2/0)) /8)/tr
Compute confidence interval width
end for Set S, E S,-1 Pir + Cirr maxjes.-1 Pgr Cjrr
Filter "bad" arms
end while return i", the only element in S, We assume i* is unique (this assumption is easily relaxed)
The algorithm proceeds in rounds, maintaining a set arms S, that are still in contention for being the best arm. At each round the algorithm pulls each surviving arm such that can construct high probability confidence interval of width c./2 around the mean of each arm. Then, any arms whose empirical mean plus confidence interval is less than the maximum empirical mean minus confidence interval are eliminated. If this the case, then that arms mean is with high probability less than the maximum mean, and is eliminated from contention. This preserves the best arm with high probability, and the algorithm terminates when only one arm remains (this best arm).
A.2 Sub-gaussian random variables
Following the exposition of [44], for strictly increasing convex function R4 + R. where (0) = 0, the v--Orlccz norm of a random variable is defined as:
X ll inf {t > O ( (X-I) < 21}.
|| is infinite if the expectation IE [bb [x)]] does not exist for any finite The sub-Gaussian parameter of random variable is defined as the 1/2 Orliizz norm, where v/2(u) e"² 1. Standard results (Hoeffding's lemma) provide that for a random variable such that a S X <
almost surely, that || || < (b-a)²
13

---

sp22 or wiitest
meds: www.gggggggaaaaaa....
hww.g Uww.a--------... en sic
11 n on kite
tur km
5.07
robk
Eestare
3.07
bur am
waki
6.065
0.000
0.000
3.000
0.155
2302
0.335
0.022
0.00
0.40
0.005
0.33
8.00
-203
6003
-200
-1re
0.300
-200
391
0,000
-1508 1001
100
520
1900 1969
Stedbrk
200g
nert am
veniamm
werd ere
testare
aod
6005
0.000
5.000
0.315
0.03
0.00-
3.800
8642
E.0L
-369
106
4mhh pul VOLUE
Am pall value
Are pal
Are (mmber
Figure 2: Distribution of arm pulls for the best, middle, and worst arms for a random query. The Gaussian fit plotted for each arm is computed using the empirical variance of the given arm, and can be seen to closely match the empirical distribution, indicating that the arm pull distributions are well approximated by a Gaussian. This merits the assumption of sub-Gaussianity Arm pulls represent importance weighted samples based on the magnitude of the query vector, so for an arm i and sample j. the value is A., sgn(x,, d. is computed across all samples for arm, and windows truncated at the lowest and highest ends of the 30 ranges across arms for viewing clarity.
Hoeffding's inequality provides a useful concentration bound, where for X with ||X|| 00² S 0²,
t?
P P(X - - 2 t) I exp
20
A.2.1 Sub-gaussianity in practice
The assumption of sub- uub-Gausiaanity is the only assumption that we make in this paper. is one of the weakest assumptions possible (does not assume that the arms are Bernoulli or Gaussian), and is a common assumption in the multi-armed bandit and adaptive computation literature [4]. Unfortunately, without such an assumption, no nontrivial results are possible; consider the case where we do not have preprocessing access to A, the vector is all 1s, and A is all is except for randomly selected entry which has value 2. In this case, any algorithm for PAC computation of softmax (Ax)) with 8 = 1 1/n (even just identification of the largest entry of As) requires M(nd) samples. More practically though, these vectors are the result of a machine learning pipeline, and not of adversarial construction. As shown by our simulations, this worst case scenario never occurs in practice, and arm pulls are generally well approximated by a Gaussian (see Figure 2). Additionally, note that for any fixed problem instance, all arm pulls are bounded, and are thus sub- Gaussian.
A.3 Improved estimators
To provide theoretical guarantees in multi-arma bandit problems, stringent assumptions are often required, e.g. Assumption 1. This is so that we can provide high probability guarantees on the concentration of the estimator constructed as the empirical mean of the observed samples. Often, these assumptions are phrased as either that random variables corresponding to arm pulls are bounded in [0, 1] a.s., or that they are o² sub-Gaussian, with a known bound on o². Such analyses have been generalized to bounded random variables with a known bound, where the algorithm is able to
14

---

adapt to the unknown variance [33, 35]. In this work, as is often done to make multi-armee bandit algorithms more performant [4], we instead directly use the empirical variance of an estimator as its aub-Gaussian parameter. In the Gaussian case, random variable's sub-Gaussian parameter o² matches its variance: as our arm pulls are constructed as the sum of many (d) terms, which can be thought of as weakly dependent, our arm pulls can be thought of as Gaussian random variables with variance o². Since in practice, we do not have a good bound on the magnitude of these arm pulls, we directly use Hoeffding's s concentration inequalities [44] with the empirical variance o², opposed to an empirical Bernstein type concentration inequality [35].
A.3.1 Randomized Hadamard Transform
We discuss applying a randomized Hadamard transform to reduce the uub-Gaussian parameter of our estimators. Define the rotation matrix R = Va DH, where D is diagonal matrix with diagonal entries equiprobably +1. H is Hadamard matrix (d must be power of 2, 0 padded if necessary). A and R to æ yields arms with better sub- Gaussian
Then, we have that applying the transform parameters. Concretely, define Z = AR and y R x. Analyzing y first, we have that entry of is too large, as:
P([lv]oe 2 () S Σ
IP [. 21)
(6)
1=1
d
= nlP
Σ
E
1/2/1/1 21
j=1
212
< 2n exp
(8)
4 4 d
j=1
= exp
(9)
|
The first inequality is a union bound over the n points and plugs in for y R σ. The second equality uses the fact that the Hadamard matrix multiplied by the random diagonal #1 matrix D makes RG i..dd +1, and we use E, to denote these Rademacher i..dd #1 random variables. Next we use Hoeffding's inequality. Finally, we simplify.
Thus, with probability at least 5. // < |2
2 1og(2n/6)
similar argument can be made for Z AR, showing that each entry in the i-th row is upper bounded by || |Alllay (2log(2nd/6) holding simultaneously for all í,
d
Since Ar = Zy. we can use bandits to approximate Zy instead of Ar. With the above analysis, we have a bound on the maximum entry of Z.ryJ. giving us a bound on the sub-Gaussian parameter of each arm with probability 2
2 llog4nnd00) 8)
max Ziz0yl = max |A|l|||=|||
(10)
d
Whereas before:
max AAyæ,1 = max :
||A ||21100
(11)
In practice, did not observe that this transformation yielded much improved variance, as opposed to simply importance sampling. Thus, we do not utilize this in our main algorithm.
A.3.2 Importance sampling
Instead of naively sampling each coordinate uniformly, we can construct improved estimators using importance sampling. Concretely, consider the unbiased estimator Z where P(Z = 1 A, = P
p;
15

---

for some probability distribution {pj} where Pj > 0 for all j and E = Now, we are left with the design choice of how to construct {pý}; naively, this is 1/d. Unpacking the variance of our estimator, we see:
Var(Z) = Var(Z)= =- - n²
(12)
Consider the case where Aj Lid. Q for some distribution Q with E[Q" X² (an empirically not unreasonable assumption), where { are fixed constants. We assume is known, but A isn't. Then, we simplify the expected variance with respect this randomness A as
E [aa 2)] = E
Σ
(Azzz)²
- 12
Pj
2
X2
12
(13)
Pj
Simplifying this for nnive where Pj 1
7. and for Zap where P. x |æyll we have
[[Var(Znaiee]] dI |=|==]] n²²
(14)
E|Var Zopi)] ||=|1 11²
(15)
In the case where we are using the same matrix over many vectors x (as is often case e.g. in make leverage scores function of
LLM inference, multinomial logistic regression), precomputed quantities based off of A, not just æ. In this case, makes sense to consider forms of leverage sampling: in this work consider taking Pj & [zy Σ, Asyl.
A.3.3 Sampling without replacement using Gumbel trick
Weighted sampling with replacement becomes wasteful in larger sample size regimes, for which the same high-weight elements sampled repeatedly. It therefore becomes desirable remove sampled elements from consideration after they are sampled, re-weighting the remaining elements accordingly. We could naively repeat this iterative process of sampling, removing, and re-weighting our elements until we ended up with a sample of the desired size, say k. However, this process is sequential and quite slow. Fortunately, as noted in [43], sampling with replacement according of weights is equivalent to perturbing the logits λ1... of our desired sample weights with draws from the i..d standard Gumbel distribution and taking the elements with the top-k perturbed logits as our sample, as detailed in Algorithm 6 and 7. This process is easily batched and is much faster as result,
Further, taking the (k: 1)-th perturbed logit j in our sample solely dependent on whether or not its perturbed logit value exceeded this threshold.
empirical threshold 7, the inclusion of element
This derivation is detailed in [27], and gives us the following expression (Gumbel CDF) for the inclusion probability of element j in set S of size k drawn without replacement according to weights
w:
J (j S)
(16)
J = (ep(-enp exp(-exple - -))
(17)
(jES)
(18)
T
These empirical estimates of the marginal probabilities of selection for each column allow to generate a sequence of estimators for each arms' mean, with improved variance discussed in the next section.
A.3.4 Variance estimation for Gumbel Samples
The Gumbel sampling trick used in A.3.3 with the fixed empirical threshold also provides a different lens on sampling process. Namely, we can compute these closed form inclusion probabilities
16

---

and by setting this empirical threshold, we may treat the inclusion of separate elements as independent. Given these observations (sampled S), an unbiased estimator for the variance of f.. constructed as the importance sampling weighted mean of the observations, can be computed from [18] as:
sample elements without replacement from [d] according to π
(19)
N4 1
1
(20)
5
J
JES
E [ūs]]
μ4
(21)
Var(A;; Σ (A
njj
(Ar. AAKO,OO))
Tj. my
(22)
T
jES
j
1- fjj
Var(p):) Σ (Ai
(23)
#2
JES
Following the analysis of [18] and assuming that our threshold T fixed, we may conclude that the estimate of the variance of üú the quantity Var T), is unbiased estimate of the true variance Var(ū:) Further, in all datasets we analyzed, the entries A were generally symmetric, zero-mean and not correlated with the corresponding entries of æt seen in Figure 3(a). Since and the values of æ are selected solely based on the values of æ, these two further assumptions make the second TT
order term As,.Aaa
Tijkk
zero expectation.
Thus, for simplicity, we discard this latter summation and treat Var(ū:) an unbiased estimate of the variance throughout our implementation, which can be computed linear instead of quadratic time (updated in constant vs linear time). In practice, as observed in Figure 4, this variance estimator
Var(A) (solid green line) provides far better estimate of the Gumbel sampler's 's true variance than other methods.
Nistgrrmm weights values
Distribution
15
10
232
72
-10
-15
8.000
Jacund 0100 Jom
(a) Sampled entries of A and x
(b) Sampled values the second order term
Figure 3: (a) Sampled entries of A and the corresponding entries of for Mistral on the Wikitext dataset. The values of A are symmetrical about and not correlated with æ. (b) Sampled values of the second order term (AijA,,
T5.A TT
B Proofs
begin by proving a standard best arm identification result, for a slightly modified version of the round-based algorithm from [16].
17

---

104
102
E2er Samole
20"3
NSE
10
500
1000
1500
2000
2509
3000
3500
4209
Sample S2e
Figure 4: Variance estimates vs. empirical mean squared error. This demonstrates the dramatic improvement afforded by improved estimators and tight confidence intervals. [imp] indicates impor- tance sampling, [wr] with replacement, [wor] without replacement, [fpe-sparse] includes the finite population correction factor.
Lemma (Best-arm identification). With probability at least 1 1- 6, Algorithm BestArmId identifies the top softmax value correctly with probability using a number of observations at most
min 32a In (2 log2 (4/4,)
A
Proof. Following the proof from [16]. since an arm's 's mean is exactly computed after pulls, we have that the best arm will be identified with probability at least - 6 requiring for arm i number of pulls at most
min (200 In dn log² (4/D.))
(24)
A?
Summing over arms yields the desired result, noting that the best arm is pulled at most as many times it takes to eliminate every other arm (i.e. the second best arm, A A2).
Note that A, S C then the second term (d) will be selected, for a sufficiently small absolute Va
constant C. Thus,
n4 < Co log (" log(d)
(25)
A?
Hence, the total sample complexity this success event is most
Tle
320 In
4n log2 (4/ A/).)) (
min
A?
0 log log(d))
min
We additionally require a lemma for estimating the mean of the best-arm in PAC sense. Lemma (Exponential best arm estimation). Sampling an arm using
32a" 32 log(2/5)
2
samples guarantees that Gir estimates e max, pr multiplicative accuracy e, with probability at least 1 - 8.
18

---

Proof. We estimate the mean of arm & after T draws using the plug in estimator Ay. For simplicity, assume I. where in the end we scale the sample complexity by S Sub- Gaussiann concentration provides that with probability at least 1 - 6,
P 144-1 S
= c/4
which in turn implies
Ar s
Me Mu c/4 S log(1 +c)
for 0 < 1. Then, exponentiating both sides yields
- 1 e 11+ = (1 - c)e S c S (1 + c)ett.
Scaling the number of samples by 8 yields the desired result.
With these two steps in place, we are now ready to tackle the novel technical challenge of this work; estimating the normalization constant the softmax. denote the softmax normalization constant as (w) (p( E: e
Proposition 2 (Softmax normalization estimation: restatement of Proposition 1), Under Assumption 1, Algorithm will, with probability at least 1 - 6, estimate fz(μ) Σn Bus to multiplicative accuracy of e, using a number of samples at most
T = =nnT + T1
-34,20 log(6n/6)) 91 32 log(Gn) (6n/6) (E 71)2
c, a;
16330 log(12/0)
2
Proof. To prove Proposition 1, we want to show that with high probability, we can upper and lower bound our plugin estimator
(1 - ()(m) S f(A) s (1 (1)))
giving us our desired multiplicative error bound, We construct several success events that collectively guarantee our bound holds, and that occur with high probability.
E₁ is the event where our estimated optimal sampling frequencies are not too far from the unknown optimal frequencies, i.e.
à Z a2/2,
BC, 1, 1,...,,. n0
On this event, we will sample arms sufficiently in the second round.
Eg is the event where all estimators ái; are within their 2-sided confidence intervals in stage
2, i.e.
li= - ull S 20 (00g(112)//// i == I,...,
On this event, bound the error in the exponentiated estimator.
Eg is the event where the first and second order errors are small, i.e.
E f(0) Σ Bu B(m - (.)) s 2 (H)
and
E B iua) (()
These two terms arise from bounding f()).
19

---

We now show that if E1 and Eg and Eg all occur, then our desired bound holds.
Lemma 3 shows that when Fs2 holds,
f(A) - Σ Shu (1 3(m (a)) + (1: -x:)9)
= [s(a) eauu -a,) Σ en 32 1n,)2
(2)
Note that if E3 holds well, the expression is upper bounded by (1 c) since each of (1) and
(2) is bounded above by (()))
All that remains is to show the lower bound also holds. We use the global inequality 1 1 S ex to lower bound
f(1)) =
Σ gh
Σ Spct
=
2 f(a) + Σ
Bun B(A - u.).
Again, under F2, (1) is lower bounded by -e/2f(n), which implies our desired lower bound (n) S f(i). Thus, all three events holding guarantees desired bound holds.
Now that we know our bound holds on the joint success event, all that remains to show holds with sufficiently high probability. To do we invoke our lemmas which characterize the probability of each event.
The probability of E, Ez, and E3 all holding is P(E, P(E,) P(E,) P(A
Lemma says that if each arm is sampled Tn times. P(E) 2 1 - 88/3. Lemma says P(E, 2 1 5/3. Lemma 5 and Lemma 6 together show (( E2 Fn) 2 1 6/3. Thus, taken together have
P(E
P(E))((EJEE)((EEEELE)
(1 8/3)/
1 1-6, 6
desired, and done.
The lemmas and their proofs follow.
Lemma 3. For n, To 1752 log( 6n for all i, we have that the event E,,
0e 2 04/2, BC < 1, == n
occurs with probability least 1 - 8/3.
Proof. By a standard Chernoff bound, with o² bound on the lub-Gaussian parameter all arms, have that with probability at least - 8/3 that for all i 1,...,.
(()--C))
à Σ, (i, -,)
E, eBns
0ig
S4g where a; Σ, suy' and C
20 lug(5/no) is the Chernoff confidence interval width constructed To
such that C, (2 and 3C, 1. To simplify constants, use that 8/ In²(2) 17. D
20

---

Lemma 4. Sampling as n4 2 To guarantees that conditioned on E1, Eg occurs with probability at least 1 8/3 and that Ez. f(β) I f()) Σ, 3(4 (.) E; (p, p4:)².
Proof. Suppose we sample each arm 74 times. Note that by a Chernoff bound on each arm,
A - pal - v 2o² (oog(6//)))))
holds on each arm independently with probability at least 1 - 8/3n, so all arms are within the two-sided bound with probability least 1 8/3.
We upper bound the plugin estimator f(t).
f(i) = Σ ei:
(26)
Σ BuL+
(27)
< Σ
em
Spi (1 B(m - (s)) + 8 (.))))
(28)
where in (28) the upper bound e < + x² for I 1.79 the event Ez, since Eg, (14= (4:) S 1/3. This is because n; 2 To, and To samples already guarantees this.
Lemma 5 (First-order error concentration). On the event Eg E1, the first order error
G
Σ
Sn B(p (t:)
satisfies
P
G
Σ
cpps
s 5/33.
(29)
Proof. First, defining E, as the failure event Gl 2 E, Lin note that
P(E,) P(Ey)E,) + P E E2 P(
implies
P(E, E2) = P(Ef) P(E2)
22P(E)).
where we use that since < 1, P(E)) 2 1/2. Thus, it suffices to show that P(E,) S 8/6
G is sum of independent sub- Gaussian random variables, each scaled by constant. Thus, we have the two-sided tail bound that with probability at least 1 - 5/6
Σ Bus B(u (4)) I 22B² logBBoo(( log 12/5)
with the sum having sub-Gaussian parameter
Ba Σ e 2pp g² 0²02/m ///
21

---

Plugging in our value of B² we find
Σ Bun Blu ai)
21og(12/6)
Σ
280 13000//2//////
1/2
< / (12((1
Σ
20p4
1
2 E cnns
1/2
< B0 / (12g(1
2
BBn
Σ
Sp.
T
23 v (x)
where in (30) use that E1, à; 00 and n. ox T22 Thus
T T2 16,3 0²c (12(12
(30)
sufficient to yield the desired multiplicative error of c/2 with probability at least 8/3.
1v20²² 3 log(6n/5)r Σj=17 To
Lemma 6. If arm i is pulled at least
times, then on the success events
'2, o5
E1, E2 where the confidence intervals hold,
S (ny - 12y) S
Buy
(31)
and Algorithm 2 will require a number of arm pulls at most
91a² log (g((6)))) (Ei=1 76)²
(32)
Σ1 ag
Proof. Bounding the second order error, we utilize the fact that we have sampled proportional to f/ T)
On the event Eg (where n, 2 7) the second order error can be bounded as
Σ 32 ou: (#:= - -2;)²²
(33)
1=1
12
20 32 Bis (6g6)///)/
(34)
i=1
n
< 20 82 (6n/8)
(35)
TT
2 202 (og(6m/0)
eBp
(36)
T
Vi
i=1
2
2020233 322 log (6n/6)))
Bu1/2
(37)
T
We want this second order error to be at most E, Mi, and so require T to satisfy the inequality below
4 32 log(6n/6) /8) (Ei e
T 2
Σ, eme.
4/20² 32 log (6n 16) (Ei 2
(38)
E, aj
22

---

Algorithmically this is not a valid T to use. since it depends on the unknown f4. However, since á,, f are close to their true values on the good event Ej, we can use these estimates. Thus, take as our budget T for this second order error:
16v
20 16
32 log 6n/6) /6) (E==1 1
(39)
E, á;
which on the event Fr is larger than (38). This random quantity, so to analyze the requisite sample complexity, we use the fact than E1, we bound the quantity in (39)
64 64v2o B² log (6n/0) (E=11 7
T
(40)
Σ,, 0i
Note 64,2 < 91 so we use this simple constant in the statement of the lemma.
This can be directly compared to the case where we only sample according to as, not Yo which would yield sample complexity of
T 2 8no log(6n/8))8²
(41)
samples. Note that since EE
7:)²
<> E; by Cauchy-Schward this is always an improvement
(up to absolute constants) up to a factor of n.
Now we can provide proof of the main theorem.
Proof of Theorem I. We utilize the adaptive approximation subroutine with error probability 8' 5/3 and error e' = e/4 from Lemma 6. With 8' probability, it requires a number of samples at most
T S 343 0 log (18n/6)))
3630 3 log 18 18m/8) (Et11))
Σ, a,
256,320² log(36/8)
Best arm identification is called with error probability complexity
= 5633. From Lemma this requires sample
320 In (1) log2 (4/4.) .d)
min
(42)
We then estimate the best arms mean using Lemma 2 to accuracy = E/4 with error probability 8 = 5/3. This requires number of samples at most
512 82 log(6/8)
(43)
By union bound, all these algorithms succeed with probability at least 1 - 8. Analyzing the multiplicative error, we have that [1/(1 + + e')] [ [1- 2s, + 2a] for 0 < S 1/2, and that (1 + (1) + 62) E (1 + (61 + S2 + 8162). Thus, these success events, the numerator is approximated to accuracy e/4, and the denominator to accuracy e/4. The denominator error converts multiplicative error of 2/ the numerator, which combines to yield an error of 3 / as < .5. This allows us to simplify the total sample complexity as
512b² 82 log(6/0)
256 (36(3)/))
T <
2 76802 2,2 log((6/0))
22
g2
23

---

Thus
T $34830 0²
logllon//000) (18n/0)n
320² In
log2 (4/4.)
min
36302 log((88) 18m/88 (E 7t)
6 <,, aj
7680 8² log(36/0)c²
2
n
log
(g)
Σ 7
SC
3²0²
n log
n
+
Σ
min
log
d
+
/=1
log (1()
A?
Ej (tj
1=1
Algorithmically, the second stage, arm needs to be pulled a number times
n4 =17 log(6n/6)
16.20 (og(6n/8) (2 96) 3/3
E,ar ài
Σj 93
+ 1633 loo (22/8) 2 ài
Σ, ājj
<173 log(6n/0)
64v 6Av2022 log( (0/g6)/))
(E, 75) 78
+
c E,
16,32 02 loglll1(2111)) (1/
2² Σ, aj
where we upper bound on the success event Eg.
Combining this with the initial To pulls, the pulls from best arm identification, and the pulls from estimating the value of the best arm, have that the overall sample complexity for arm is upper bounded as:
n. = 34,320 log(6n/8)
=
64 32 (og(6n/8)
(E, 15)
E/ aj
16 log( 12/8)os
2²
E, aj
32a In 12n 8
log? (4/ A.))
512 3 log(6/8)
1 1}
And the total algorithmic sample complexity the success event is
<
min (7, d).
(44)

---

B.1 Interpreting the results
We work to provide a simplified (looser) bound on the sample complexity when the minimum gap is A. The worst case sample complexity in this case is when the best arm has mean P1, and all the rest have mean μl A.
This allows us to simplify the overall sample complexity as
T T
o²
3n
nog
(
12 log
+
A?
log
B
log
"Ejaj 0j
1=1
Σ 7
+ B logí (og(1/6)
82
(45)
i=1
SC
8
log
ni (
(1
1 1 n2 1
)
)) + 00000==0222 log(1/0) + n log
n log d
A
1 BD
where we use the fact that for n >
Σ 75
+ (n 1)e RA/22)
Σjag aj
1+ (n 1)e-88 13
1 + (n 1)2e - 3D
2
14 (n l)e-9A
+ 1ite-9a 99
0
1 + ne ne-
Evaluating our sample complexity in (45), when A and the sample complexity can be more simply bounded as.
2 log n, this term is bounded by constant,
log log
T T 80 log
(5)
n (n + + + :000= log(110) + no²
log n
S CC80 (log
n (n +)=+++++= 20gg(100))
(46)
In the last line we require the condition that
log
log d
3 log
n1
(47)
log n
iee. is not doubly exponential n and 3 not too small. Assumption 2 (large gap. moderate 8, moderate d). We assume that (47) holds, and that A log
Under the conditions of Assumption 2, the sample complexity in Theorem I can be simplified in Corollary I.
B.2 Asymptotic optimality of sampling frequencies
Following the approach of [8], show the asymptotic optimality of our sampling frequencies (sampling proportional to a for minimizing the first order error, and 7/4 for minimizing our bound on the second order error).
B.2.1 First order frequencies os;
Considering plug-in estimator fo; analyzing the first order taylor expansion of error, we have that
f - f(m) (n)) (A
n1) + (IlÆ (1パí-111))
Thus in the high accuracy regime (e - 0), we can consider only the first order term. Assuming Gaussian noise in our arm pulls (an identical result holds for sub-Gaussian noise), the first order
25

---

error can be bounded as (assuming we use T pulls, and sample arm i, p;T times for probability distribution p:
Vf(H) (p - 1)) ~ ~N N (0, <
Tp:
Optimizing over probability distributions using Sion's 's minimax theorem and strong duality as in [8] gives us that
20p4
a- argmin
82
Tp:
2pp
+ A EEii-
argmin max
Tp;
utilizing strong duality, we have that
0
200 + +x=
a'
x eau.
(48)
Pp:
Note that in the limit as e 0, multiplicative and additive error objectives are equivalent.
B.2.2 Second order sampling frequencies %
With the first order error term hand, the second order error term left to be analyzed:
Σ e (Mi pa:)²
(49)
We minimize bound on this given by our confidence intervals 144 -u--is A < vap. where e is some constant. We identify what sampling distribution p minimizes this second order bound, defining n4 = T.
By similar argument as for the first order analysis:
f2 8pc 1 Tp:
7*
argmin
p
Spc 1
Ep:-1
argmin max
+ X
p
Pi
utilizing strong duality, we have that
i
8r Di
-
Y
2112
(50)
Op:
Gains by sampling according to % Sampling according to gives second order error bounded
by
2
Σ
euxl2
|
(51)
as opposed to sampling according to a which gives a second order error bounded by
Σ sr. - nl|r1
(52)
By standard norm inequalities, sampling according is always at least as good, and up to factor of n improvement the case where one entry in the softmax much larger than the rest; exactly the case of interest.
26

---

B.3 Comparison with [8]
In [8], the problem of estimating real valued function additive accuracy e with probability at least 1 is studied, under the assumption that the function has L-Lipschitz gradients. For the case of softmax estimation, the gradients are not Lipschitz due to the unbounded nature of the exponential. However, if we evaluate the norm of the gradient at a point μ., we obtain
| [ = lim max c ( - vS(u))1
c-0 Sulle
lim max
Σ (e - -Be
c- Sullee
lim max c ²3²
2814 (c - -1)
e-0 Sullse
5 max e 20p.
Theorem 1 of [8] states that the number of samples required to achieve additive emror with probability least 1
V f log(1/8)
n.² L
T=
(53)
g²
where the noise variance x² is assumed to 1. Since the error in our setting multiplicative, we are interested in s' ef( = E; e9t4, Additionally, f = 3 e Thus, the number of samples required is
V [ (oo(1/))
n² m²Llog(n//)
==0
+
6²
B²
(E,
eB4 log(1/5) n² max, 2a log(n/o)
E2
EE
e844)
E, e
log(11/)
n²3t max, 2014 log(n/8)
=0
22
> eBu
This is be compared with the sample complexity of the proposed algorithm in this paper, for the specific setting of softmax normalization estimation, which is
log (1)
12 T CB2 ni nog
(
+
log
T
12
Σ i=1
Σ aj
+
x²
taking の² = 1 to compare results.
The constant term independent of linearize the exponential. The E term matches between the two settings, as asymptotically the optimal strategy is indeed to sample according to the first derivative. The term scaling with - improves dramatically on that of prior work. Note that in the case of f1) || the second order term (scaling with = can be improved to 3/, as the mean of the second order error can be removed. Thus, we can see the massive improvement afforded by more refined algorithm, tailored for the specific structure of the softmax function.
B.4 Extension to heterogeneous variances
Adapting bandit algorithms to settings with heterogeneous variances has been done in both the standard regret [3] and best arm identification [33] settings.
For best arm identification, sacrificing log factors for the sake clarity, empirical-Bernstenn--beee Bernstein-based confidence intervals [35] can be constructed where iteratively pull each arm once, try and
27

---

eliminate, and progress. Union bounding over the nd possible pulls naively upper bounds this, yielding complexity of
Σ min ((
2
A?
+ A
log
(54)
nd
This assumes that all arms are bounded in [0, 1] with variance of-
We additionally require a lemma for estimating the mean of the best-arm in PAC sense. Lemma (Exponential best arm estimation). Sampling arm i
T 32o 3 log(2/8)
22
samples guarantees that eits estimates eBs multiplicative accuracy e, with probability at least 1-
This trivially follows from the proof of Lemma 2.
For the softmax normalization estimation, we know from [8] that the optimal first order sampling frequencies are to sample arm number of times proportional to
n4 O o
(55)
However, sampling like this yields additive first order that scales E, ore which cannot be easily related to Σ; eBii, as would need to get multiplicative error bounds. Thus, instead use suboptimal target first order sampling frequencies, scaling with of. to avoid this analysis issue
n4 x 0u.
(56)
Scaling the number of pulls for each by of yields:
Proposition 3 (Softmax normalization estimation: heterogeneous variances variant of Proposition 1). Under Assumption 1, Algorithm will, with probability least 8, estimate fø(H) E, 8My to multiplicative accuracy of e, using a number of samples for arm at most
n; 34.52 log(6n/0)) (0)n 91a 82 log log(6n/5) (E 7 16,3 of log (1/66)
c.
The proof of this proposition follows similarly to Proposition 2, sampling proportional to of cancels the differing variances, essentially reducing the problem to the homogeneous setting (suboptimally).
C Implementation Details
In this appendix, present the implementation details of algorithm. In Algorithm 4, provide pseudocode with greater detail about our implementation of the Adaptive Softmax algorithm. We note that Algorithm 4 contains some implementation differences from the original Algorithm 1 presented in Section 4. None of these changes materially affect the output of the algorithm; nonetheless, provide a discussion of them here enable reproducibility of our experimental results. Our results are also reproducible via 1-line script in our code submission. In the following subsections, we describe each of these implementation details.
As an important note, we consider all variables global unless stated otherwise. This is to say, calling "pull arms" updates the state of arm mean estimates fi, their variance estimates ô, the number of pulls per arm {n:). inclusion probabilities π, etc. We define the estimators based on Gumbel sampling simply
according to importance weights the set of arms A. This idea of treating sequence of estimators with confidence intervals was pioneered for the computational setting in [4], and saw further usage in [24].
C.1 Reusing Arm Pulls
In our theoretical analysis Appendix B, each phase of Algorithm handled independently. This allows to union bound the error probabilities of each phase of the algorithm. In our implementation,
28

---

Algorithm 4 Adaptive Softmax (implementation details)
1: Input: Matrix A, vector 2, temperature 3, error failure probability 8
2: Output: With probability at least - 6, the argmax coordinate i" and an estimate p:- of its
probability such that c)p-- v P;- I (1 + c)p.-
3: w < Get.mportanceeeiigiee
(A
b Algorithm
4: P,
5: Construct set of arms A from A, x, P.c
6: o² =
Initial variance to pull arms to
7: PullToVariance
>> A Agorithm
8: 1" < (A.
> Algorithm
9: P... + (A,,.)
Exact computation of f44 = μ.
10: Z Normal sarraaaiaattatttiiniiimaaa......
> Algorithm
11: Compute estimated probability as fA- e /2
12: return i" Pr-
Algorithm 5 Get GetIImportanceliiggin
Input: Matrix RRndd vector x E Rd
2: Output:
Rd vector of importance weights
3: for all = ,,...,d do
4: Compute w [r Σ D==1 [Asyl
l norm of ith column of A
5: end for
6: return
Element- multiplication of l, norms of columns of A, and [ac]
Algorithm 6 Gumbeeleermuattion
1: Input: Importance weights w
2: Output: Permutation P, cached outputs for inclusion probability calculation
Lid
3: Draw each & ~ Gumbel(0, 1)
4: L log w
Log importance weights D Log importance weights perturbed by i.i.d. Gumbel noise Compute thresholds in decreasing order Compute sorting order of thresholds (argsort) Store log importance weight and sorted thresholds for later use
5: + ++
6: h sorted(L')
7: P + (Ldering(L')
8: c c+ hL
9: return P.c
Algorithm 7 InclusionnPoobabllttes
Input: Sample size k; cached outputs c (4, h) from GumbelPermutation
2: Output: Inclusion probabilities vector
3: L, t-c
4: H = -00
> Initialize cutoff threshold
5: k then
6:
И + h
> Set cutoff threshold to kth largest perturbed log importance weight Inclusion probabilities; derived from Gumbel CDF
7: end if
8: π = 1 - -exp(- exp(I. H))
9: return
Algorithm 8 PullToVariance
Input: Set of arms (sequence of estimators) A, variance a²
2: Set Ç .1, multiplicative pull increase factor
3: while there exists an arm with o² > o² do
4: A' - {i E A o7 > o²), n, is corresponding number of pulls
5: PullArms (1 + ()n,)
6: end while
29

---

Algorithm 9 PullArms
Input: Set of arms (sequence of estimators) A, target number of pulls per arm {Ni
2: for all arms do
3: if n4 N, then
4: continue
Do not pull arm i if it has already been pulled N; times
5: end if
6: Compute h; = -;
7:
n' < pacllusiionrrlbbbliiieeee..:
8: Sample k; coordinates without replacement according to weights n'' X
9:
Update mean estimate fb,
p4 + ky.
E==1 si
N:
N;
2 o + N? E.= (AL 21:5( N?
10:
Update variance estimate o²
11: end for
Algorithm 10 EstimateArs
1: Input: arm i A, pull variance of, multiplicative error e, failure probability
2: Output: Mean estimate /4s
targetVar< B0
4: PullToVariance targetVar)
5: return 14
however, we re-use arm pulls across different parts of Algorithm 1. Intuitively, once an amm return has been observed (corresponding to scalar multiplication of an element of A with an element of r), it can be used to warm-start estimates of ji and 6² in later stages of the algorithm. In practice, we observe that the re-use of arm pulls does not affect the correctness of our algorithm and yields significant sample complexity improvements compared to cold-starting each stage of the algorithm independently.
C.2 Exact Computation of Best Arm
In Line 9 of Algorithm 4, we compute the mean of our estimated best arm i" exactly, and set the "estimate" e = Pr- This allows us to reduce the approximation fidelity required by Normalizationnssiiaati
(Algorithm 2) to instead of and saves a constant factor in sample complexity. Furthermore, this computation of μ4 efficiently computed as a vector-vector dot product.
In practice, found that Algorithm usually required over d samples for the best arm. As such, performing the computation /4 2x (reusing coordinate-wise samples from previous stages of the algorithm) did significantly increase sample complexity.
C.3 Initial Pulls (To)
(Algorithm 2), the initial number of arm pulls To depends on o². However, as discussed in Appendix A.3, this variance proxy is often unknown a priori. As such, we set To = IO
In
d. We observe that this choice of To works well experiments, across multiple datasets.
C.4 Tuning
We note that despite the changes made above, there still exists some looseness in practice. To remedy this, we scale our variance estimates by a constant factor to reduce the amount of pulls needed to reach target variances, and improving the gains of our algorithm as result. We generate these constant factors for each dataset/model by tuning on a separate "training" group of queries. Tuning is performed, generally, via bisection to discover failure probability parameter 8. This bisection is performed in geometric space and terminates when the log difference between the low and high end of our interval is within 10-2, The range of factors consider is [10-6, 1]. We first tune the constant factor independently for the variance estimates
minimal factor which still satisfies our provided
30

---

used in bandits such that the algorithm successfully identifies the best arm with rate of at least 1- on our training set. Next, we tune the constant factor for the variance estimate used in log norm estimation such that the entire SFTM algorithm succeeds with rate of at least 1 - å on our training set.
C.5 Wall-clock improvement
The focus of this paper is develop the first provably adaptive softmax algorithm with PAC guarantees, highlighting its dramatically improved sample complexity across a wide variety of models and tasks. The eventual goal of this method is to enable wall-clock improvements in hardware implementations. These next steps of converting our provably and empirically performant method into a hardware optimized wall-clock efficient algorithm is an exciting direction future work, which we detail below. In most modern-day transformer architectures, memory I/O serves as the primary bottleneck [20]. AdaptiveSoftmax already presents an opportunity to significantly scale down the number of entries of the matrix that must be loaded at inference time, and, in the future if memory remains the bottleneck improve model bandwidth by a similar factor. This objective appears in reach, since we have designed the components of AdaptiveSoftmax to be amenable to tiling and parallelization. Most notably, our implementation of AdaptiveSoftmax uses the same column to generate an decision are two-fold. First, takes advantage of the locality of entries the same column to load samples faster, and, second, it removes intra-column correlation, which can yield theoretically improved performance [7]. Adjacent column samples can also be combined by simply summing their respective importance weights admitting simple tiling of our matrix that could easily be sized particularly to fit individual tiles into SRAM on a GPU along with a copy the vector and the current mean/variance estimates for each arm. Then, based on the arm estimates as we do currently. The successive elimination bandits algorithm utilized by AdaptiveSoftmax is also, by choice, quite easily parallelizable. We may also store two copies of our matrix one with wider tiles and one with taller tiles - to take advantage of this structure at all stages of the AdaptiveSoftmax algorithm: both when larger number of samples is necessary for fewer arms, in later stages of adaptive estimation, and when a smaller number of samples is necessary for many arms, in earlier stages of the adaptive estimation. This said, we observe in our experiments that the bulk of compute is invested in our early samples of many arms. Just using basic parallelization to speed up this step could result in the desired speed improvements.
weitghtee sample for each active arm. The reasons for this implementation
dynamically load these tiles into SRAM
C.6 Robustness to parameters
The user-desired parameters can take a wide range. We kept 30% constant across all simulations because we observed varying did not result significant changes to the performance of adaSoftmax. Further, we suspect that for most users, 5 will fall in the range considered: 90- 99%. However, to assuage any concerns and verify our assertion that adaSoftmax is not sensitive to choice of E or ó, we include here adaSoftmax with much wider range of parameters on the MNIST dataset in Table 3.
0.001
= 0.01
8 1.0
à 0.01
100%, 5.18x
6.63x 100%,
99.38%, 8.13x
99.38%, 8.14x
8 0.05
99.75%, 6.64x
99.13%, 8.46x 95.38%, 8.80x
93.50%, 8.81x
8 8= 0.2
91.25%, 7.19x
89.25%, 8.90x 88.00%, 9.29x
83.375%, 9.25x
Table 3: Success rate (%) and FLOP gains (x) for adaSoftmax with varied o and on the MNIST
dataset, showing the improved performance across wide range of parameters, and that raising past 0.1 causes minimal difference in performance.
D Additional extensions and comments
D.1 Effect of temperature
Temperature is treated as fixed constant (fixed parameter for the problem at hand, not tunable by the algorithm). This is because tuning the temperature fundamentally changes the problem. With higher temperatures, the only arms that matter are the best and second best arms, and adaptivity
31

---

is extremely helpful. At low temperatures, the output will be essentially the uniform distribution, and the computation is trivial and adaptivity unhelpful. With respect other parameters, the error probability and FL.OP gains of adaSoftmax are insensitive to changes in and vary most with the choice of 6. We demonstrate this trend on the MNIST dataset Table
D.2 Application Nucleus Sampling
The analysis posed in this paper focuses on identifying the largest entry in the softmax output and estimating its associated probability. As discussed above, this naturally extends to identifying the largest elements in the output vector by replacing the bandit best arm identification algorithm with any top-k identification algorithm [39]. However, LLM inference, the goal is often draw sample from the softmax output distribution via nucleus sampling [17]. Nucleus sampling avoids specifying k directly; instead, provides cumulative probability p and requires the identification of the top elements such that & is the smallest value such that the sum of the probabilities of the top k elements greater than p. The next token is then sampled from the renormalization probability distribution on these k elements. Our adaptive sampling algorithm naturally applies to the nucleus sampling setting. AdaptiveSof tmax can maintain a predicted set arms such that the sum of the arm probabilities is greater than based pessimistic S, sampling both the arm with the lowest mean minus LCB (in an attempt to verify the boundary), as well as the arm with the widest confidence interval (in order to better estimate p), and b) sampling the top arm in [n] S. to see it belongs in S. For simplicity and concreteness, in this work we focus on identifying and estimating the probability of the top-1 1 element, but this is exciting direction of future work.
mean estimates. Then, iteratively sample: a) arms in
32

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
The abstract and introduction clearly state the theoretical and practical contri- butions of the paper, highlighting its theoretical guarantees, sample complexity reduction, and potential for wall-clock improvements.
Justification:
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
We discuss assumptions on known uub-Gaussian parameter bounds Assump- tion I and computational difficulties in realizing sample complexity gains wall-clock gains.
Justification:
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [Yes]
Justification: All formal claims are stated with their requisite assumptions in the main text, with proof sketches. Detailed theoretical proof's are provided or cited in the appendix.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes)
Justification: We provide a one-line reproducibility script to reproduce the results in the paper.
5. Open access to data and code
Question: Does the paper provide open access the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: Our codebase publicly available on github, and is reproducible via a 1-line reproducibility script: attps://giis ctt. comm/hruuGroup// wwwaaaiapsasooooe tmax.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Answer: [Yes]
Justification: This is described in the paper
7. Experiment Statistical Significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: Error bars are provided for all plots.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
33

---

Answer: [Yes]
Justification: This is described in the paper.
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics https: : htt///nerrppppp hwwprrpssccppubi///tttiiieeee hccppbblcc/thhiccsuuee1l111111
Answer: [Yes]
Justification: The research in this paper conforms with the stated Code of Ethics.
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: We discuss broader impacts in the final section of the paper. There minimal societal effects, and potential energy and environmental savings from more efficient computation.
11. Safeguards
Question:
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Answer: [NA]
Justification: The algorithms proposed this paper do not have high risk for misuse.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: All resources are publicly available.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: We provide a publicly available github bectps///gitin com/ThrunGroup/ adaptiveSof tmax that is well documented and has 1-line reproducibility script.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: This paper does not use or perform any research with human subjects.
15. Institutional Review Board (IRB) Approvals Equivalent for Research with Human
Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: This paper does not have any human subjects, and does not require IRB approvals.
34


Reviews:
[{"review_id": "r2U57Ev77L", "reviewer": "Reviewer_qEfP", "summary": "This paper introduces an algorithm named AdaptiveSoftmax, designed to efficiently compute the top-k softmax values rather than the full softmax computation. This algorithm is particularly useful in high-dimensional settings where the computation of the full softmax function can be prohibitively expensive. The authors provide theoretical guarantees for the efficiency and accuracy of AdaptiveSoftmax and demonstrate its effectiveness through empirical results on real and synthetic datasets.", "questions": "a)See weakness", "limitations": "8.The limitations is discussed in the paper.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "a)The paper presents a novel approach to approximating the softmax function using adaptive sampling.\nb)The authors provide probably approximately correct (PAC) guarantees for the performance of the AdaptiveSoftmax algorithm, which strengthens the credibility of their approach.\nc)The paper includes extensive empirical results demonstrating the algorithm's efficiency on both synthetic and real-world datasets, including significant reductions in sample complexity compared to brute-force softmax computation.", "weaknesses": "a)The algorithm relies on assumptions about the data distribution and variance proxies, which may not hold in all practical scenarios.\nb)The adaptive nature of the algorithm introduces implementation complexity, particularly in balancing computational resources and ensuring efficient sampling.\nc)While the algorithm is theoretically sound, its practical benefits are most pronounced in scenarios with very high-dimensional data. In lower-dimensional settings, the gains may be less significant.\nd)It is better to give detailed explanation of each step in Algorithms. What is EstimateArm in Algorithm 1 step 8?\ne)How sensitive is the algorithm's performance to the choice of parameters such as the temperature parameter?\nf)Besides comparing with the full softmax computation, is there a detailed comparison with other approximate softmax algorithms?\n\nTypo:\nLine 5: “we present present an …” should be “we present an …\nLine 135: “Our objective them becomes …” should be “Our objective becomes …”", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer qEfP for their insightful and detailed feedback.\nWe provide a point-by-point response to Reviewer qEfP’s comments below.\n\na) **Sub-gaussian assumption:** This assumption is minimally restrictive, and is borne out in practice.\nWe provide an in depth discussion in main response point 1, which we have added to the Appendix.\n\nb) **Computational efficiency of adaptivity:** our algorithm provably obtains sample complexity improvements, and as discussed in main rebuttal point 3, there are several exciting next steps to a final hardware optimized implementation that realizes these wall clock gains.\nOn the theoretical side, there has been a surge of work over the last decade focused on developing multi-armed bandit algorithms with minimal rounds of adaptivity that still retain the same theoretical guarantees.\nAs we show in this work, normalization estimation can be accomplished in 2 rounds of sampling, best arm identification in $\\log(1/\\epsilon)$ rounds, and best arm mean estimation in 1 final round.\nThe round complexity of PAC best arm identification can be further improved, at the cost of worse dependence on the gaps ([13], Hillel et al.).\nMain rebuttal point 3 provides additional discussion, and many of these notes can be seen in our publicly available and easily reproducible code base (see attached files, will be made public on GitHub after the paper is published).\n\nc) **Dimensionality:** Our algorithm indeed displays its largest gains in the high dimensional regime, which is the natural regime future models will be focused on.\nAs shown by our numerical experiments though, these gains are still substantial for moderate $d$.\nFor example, considering GPT2 ($d=768$), one of the smaller language models in use, our approach still yields over 7x improvement in sample complexity.\n\nd) **Algorithm clarity:** Thank you for the suggestion.\nWe were space-limited with the initial submission and moved the algorithm descriptions to the Appendix.\nWith the added page for the camera ready version, we will add these algorithmic descriptions back into the main text.\nThe function EstimateArm simply pulls the input arm to accuracy epsilon with probability at least $1-\\delta$ (as in Lemma 2): we will clarify this and package the result as an algorithm in the final version.\n\ne) **Temperature:** Temperature is treated as a fixed constant (fixed parameter for the problem at hand, not tunable by the algorithm).\nThis is because tuning the temperature fundamentally changes the problem.\nWith higher temperatures, the only arms that matter are the best and second best arms, and so adaptivity is extremely helpful.\nAt low temperatures, the output will be essentially the uniform distribution, and the computation is trivial and adaptivity unhelpful.\nWe will add a brief discussion clarifying this to the Appendix.\nWith respect to other parameters, the error probability and FLOP gains of AdaptiveSoftmax are insensitive to changes in $\\epsilon$ and vary most with the choice of $\\delta$.\nWe demonstrate this trend on the MNIST dataset in Rebuttal Table 1.\n\nf) **Comparison:** As we discuss in overall point 2, there are no other approximate softmax algorithms that can provide PAC guarantees aside from the baseline of exact computation, which provides our point of comparison. We provide a qualitative comparison to a baseline without PAC guarantees suggested by another reviewer kEa9 (see **Comparison**).\n\ng) **Typo:** Thanks you for catching this, fixed.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "GppSt3VXnR", "reviewer": "Reviewer_qUNV", "summary": "The softmax function is a widely used tool, e.g., as an activation in the final layer of a classifier. Hence, cutting down on its computational costs can have a significant impact across the AIML field. This paper aims at this by introducing an adaptive variant of computing the softmax for the top $k$ values by estimating the normalization and the index of the likely highest value in the input by employing a multi-armed bandit setting. Furthermore, theoretical guarantees on the accuracy of the output of the adaptive algorithm are provided and can be controlled by an additional parameter $\\delta$.", "questions": "I have no questions for the authors.", "limitations": "The authors have addressed the limitations in a dedicated section.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "- Considering how ubiquitous the softmax function is in AIML methods, this paper promises a significant impact. The reported gains in sample efficiency within both synthetic and real-world datasets further underpin this point.\n- The method has been presented with good clarity and has a good level of originality. Creating more flexible computation models is an important research endeavor when tackling the growing resource demands by deep learning methods while maintaining a desired level of accuracy, which is honored in this work by allowing for fine-granular control of the trade-off between resources and accuracy by using the target error and failure probability parameters.", "weaknesses": "- Overall, a more extensive evaluation with improved clarity is highly desirable. It would be important to see the effect of changing the $\\epsilon$ parameter as well when conducting the evaluation, which is, to my understanding, set to a constant value of $30%$ throughout the experiments. Furthermore, reporting the effects of these changes could be made more comprehensive if presented in a plot with a wider range of choices for $\\delta$, rather than just the three.\n- As has been discussed in the limitations section, there can be a considerable trade-off between more involved, adaptive methods and easy-to-batch and parallelize brute-force computations. To my understanding, the experiments do not consider this trade-off, e.g., by reporting the wall-clock time next to the gains in sample efficiency.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer qUNV for their helpful feedback regarding improving the exposition of the improved algorithmic performance of our method.\nWe respond to their two main points below:\n\n1. **Varying parameters:**\nThe user-desired parameters can take a wide range.\nWe kept $\\epsilon = 0.3$ (i.e. 30\\%) constant across all simulations because we observed varying $\\epsilon$ did not result in significant changes to the performance of AdaptiveSoftmax.\nFurther, we suspect that for most users, $\\delta$ will fall in the range we considered: $0.1-0.01$ (equivalently, 90-99\\%).\nHowever, to assuage any concerns and verify our assertion that adaSoftmax is not sensitive to choice of $\\epsilon$, we now run adaSoftmax with a much wider range of parameters on the MNIST dataset in Rebuttal table 1, and will add this figure to the appendix.\n\n2. **Converting FLOP gains to wall-clock gains:**\nAs discussed in main point 2, our algorithm can be made minimally batch adaptive, to improve its parallelizability.\nWe further discuss in main point 3 the exciting line of future work towards converting the significant reduction in FLOPs achieved by AdaptiveSoftmax to improvements in wall clock speed.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "iSo94j0Ppw", "reviewer": "Reviewer_kEa9", "summary": "This paper focuses on the efficient approximation of the softmax function. The authors propose an algorithm named AdaptiveSoftmax, which aims to reduce the computational cost of the softmax function in high-dimensional environments. Inspired by the multi-armed bandit problem, the algorithm adaptively allocates computational resources to important output values, efficiently calculating the top k values of the softmax function. It also provides PAC (Probably Approximately Correct) guarantees, demonstrating its effectiveness both theoretically and experimentally.", "questions": "How was the performance of AdaptiveSoftmax evaluated? For example, are the speedup metrics in Tables 1 and 2 theoretical values or execution time comparisons? If it is the latter, what benchmarks were used? PyTorch’s torch.nn.CrossEntropyLoss utilizes log_softmax and nll_loss, which are highly optimized for GPU performance. Has AdaptiveSoftmax been compared with these in practical settings on GPUs?", "limitations": "The authors appropriately mention the limitations of the proposed algorithm, noting that it is most beneficial in high-dimensional settings. However, in practical settings involving LLMs, methods highly optimized for GPUs are usually used. Therefore, the question remains as to whether a fair comparison has been made regarding the current effectiveness of the proposed method in this setting.", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "Originality:\n\nThe AdaptiveSoftmax algorithm introduces a novel approach by adaptively allocating computational resources, addressing challenges in existing softmax approximation methods. Traditional approaches, such as using hierarchical models, reduce computational complexity but increase the number of internal classifiers, lacking accuracy guarantees. Methods leveraging label prior distributions require prior knowledge, limiting their applicability.\n\nQuality and Presentation:\n\nThe paper provides detailed theoretical analysis and experimental validation, robustly supporting its technical claims. The presentation is clear, with a well-organized narrative from the review of related work to the proposal of the algorithm, its theoretical guarantees, and experimental results.\n\nSignificance:\n\nThe paper offers a novel solution to a significant problem in computing the softmax function. The AdaptiveSoftmax algorithm not only substantially improves computational efficiency compared to existing methods but also provides unique PAC guarantees. This algorithm holds potential for significant computational efficiency improvements in machine learning models dealing with high-dimensional data, potentially impacting a wide range of applications.", "weaknesses": "In practical implementation, parameter tuning may be necessary. An ablation study on the sensitivity and justification of default hyperparameters would enhance the paper. Additionally, discussing scenarios where the proposed method may not be suitable would strengthen the paper. For instance, while the performance degradation due to approximation is shown to be limited on average, what about the worst-case scenarios, such as cases with significant data distributional shift?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer kEa9 for their in depth review: we wholeheartedly agree that this method of adaptive computation holds the potential for significant computational improvements across a wide range of applications. We respond to the specific questions below.\n\n1. **Ablation study:** The theorems related to our algorithm’s sample complexity suggest minimal dependence on $\\delta$, and slow growth with $\\epsilon$ in the moderate $\\epsilon$ regime, which is consistent with what we’ve observed experimentally. Following your suggestion, to better demonstrate and highlight this point, we have now added a table showing the sample efficiency scaling for a wide range of values of $\\epsilon$ and $\\delta$ in the Appendix (Rebuttal Table 1).\n\n2. **Limitations:** As we discuss in the Limitations section of our paper, this method attains its largest gains in the high dimensional setting. However, as can be seen from our experimental results, even in smaller models like GPT2 ($d = 768$) our model exhibits over 7x improvement in sample complexity.\n\n3. **Distribution shift:** Fortunately, our method does not suffer from these issues. It is designed to run at inference time and is truly instance-adaptive to the weights of the underlying model and the query. This is a great strength of our method, because as we discuss in main comment 2, many existing works such as the one referenced by Reviewer wft8 do in fact suffer from issues relating to distribution shift. Thank you for highlighting this point, we now add additional discussion to this end in our related works.\n\n4. **Performance evaluation:** in this work, we use a proxy for FLOPS as our metric of comparison. This works by counting how many entries of the matrix $A$ the algorithm needs to observe (essentially, how many multiplications need to be performed). As this work focuses on providing a novel softmax approximation method with provable guarantees, with the goal of minimizing the number of entries of $A$ that need to be observed, this is the relevant metric to use to see that our algorithm is yielding the gains predicted by theory. The baseline that is compared to is brute force exact computation, which computes the entire $n \\times d$ matrix vector product. Our algorithm requires many fewer samples across the wide range of settings and parameter regimes we tested.\n\n5. **Comparisons:** The suggested methods all perform exact computation, and so are captured by the naive baseline in our comparisons. These methods are highly GPU optimized, implementing considerable hardware-conscious optimization, and so while their sample complexity is much worse than adaSoftmax, their wall clock complexity is better. Given the significant sample complexity improvement afforded by our method, we believe that improved wall clock performance is imminently possible, as we detail in Future Work. We discuss these details in points 2 and 3 of our main response.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "BUjxFamL40", "reviewer": "Reviewer_wft8", "summary": "This paper introduces an efficient algorithm called AdaptiveSoftmax, designed to compute the top k outputs of the softmax function more effectively than the traditional full softmax computation. The key innovation lies in its adaptive approach, which leverages multi-armed bandit techniques to prioritize computations for the largest input elements, significantly reducing sample complexity. The paper provides PAC (Probably Approximately Correct) guarantees for AdaptiveSoftmax, demonstrating its theoretical soundness. Empirical results on both real and synthetic datasets, including the EuroSAT dataset, validate the algorithm's efficiency, showing substantial reductions in computational overhead, often achieving tenfold improvements or more. Additionally, the proposed method for estimating the softmax partition function offers potential applications beyond the current scope, such as in kernel density estimation.", "questions": "1. Please discuss if the assumptions are a constraint for the environment, does it have an impact on the effectiveness of the algorithm?\n2. It would be nice to see some quantitative comparison of the proposed method with other existing algorithms, eg., speed improvements.", "limitations": "The authors have provided an extensive explanation of the limitations in the final section of the paper.", "rating": 5, "confidence": 1, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1. The Adaptive Softmax algorithm significantly reduces computational overhead compared to the full Softmax computation. This is particularly beneficial in high-dimensional settings, making the approach practical for large-scale Machine Learning applications. The concept of focusing on the relevant top-k outputs is interesting.\n\n2. The paper provides strong theoretical foundations with PAC guarantees for the AdaptiveSoftmax algorithm. These guarantees ensure that the method is not only empirically effective but also theoretically sound, offering reliable performance bounds. There have been extensive experiments on real and synthetic data. The authors also provide results on a variety of networks, as CNNs and LLMs.", "weaknesses": "1. The method relies on the assumption of a variance proxy bound for the sub-Gaussian parameters of the constructed estimators. While the paper discusses loosening this assumption, its practical implications and the extent to which it holds in various scenarios are not thoroughly explored, potentially limiting the algorithm's applicability in more varied or less controlled environments.\n\n2. The paper does not compare quantitatively to other adaptive Softmax methods, like A\n\n[A] Joulin, Armand, Moustapha Cissé, David Grangier, and Hervé Jégou. \"Efficient softmax approximation for GPUs.\" In International conference on machine learning, pp. 1302-1310. PMLR, 2017.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer wft8 for their careful reading of our discussion of empirical and theoretical improvements, and for highlighting this relevant reference. We discuss this paper, and the assumption of sub-gaussianity below.\n\n1. **Sub-gaussian assumption:** this is a very fair question, which we discuss in detail in main response point 1. Briefly, assumptions are necessary to avoid $\\Omega(nd)$ sample complexity, and empirically our assumptions are borne out. Folklore suggests that trained weights of large models are generally normally distributed, and we have verified this for the models examined in this paper (see attached Figure 1).\n2. **Comparison:** following up on main point 2, no other methods provide theoretical guarantees (our PAC formulation in equation 4) aside from the naive baseline of exact computation, which we compare to. The method of the suggested paper \"Efficient softmax approximation for GPUs\" does not provide PAC guarantees, and only gives empirical wall clock speedups without any theoretical analysis of the number of samples needed, making it difficult to compare directly. At a high-level, the suggested method finds a learnable quantization (referred to as ``clusters'') of the matrix $A$, where the weights of the clusters are learned during training. Already, this poses a problem for inference since the vocabulary in the test set may be distributed differently than the training set for which the clusters are optimized. This falls into the concerns made by reviewer kEa9 of distributional shifts which our algorithm does not suffer from. We attempted to implement this method to compare against, but ran into some issues with their implementation. For example, note that the *forward* method (defined as a *torch.nn* module at [link](https://pytorch.org/docs/stable/generated/torch.nn.AdaptiveLogSoftmaxWithLoss.html) requires as input *target* which are the output labels that are not available during inference time. In comparison, our algorithm includes the cost of approximating these labels (i.e. the best-arm identification) and does so efficiently while providing PAC guarantees.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We would like to thank all the reviewers for their careful reading of our manuscript.
We were pleased to see that all reviewers appreciated the novel PAC guarantees provided by this work for efficient, instance adaptive softmax computation.
We have addressed all the comments and suggestions made by the reviewers, which has helped improve the quality and clarity of the paper.
We discuss some common points of concern below, and look forward to the upcoming discussion period.

1. **Assumption of sub-Gaussianity:** This is the only assumption that we make in this paper.
It is one of the weakest assumptions possible (does not assume that the arms are Bernoulli or Gaussian), and is a common assumption in the multi-armed bandit and adaptive computation literature [4].
Unfortunately, without such an assumption, no nontrivial results are possible; consider the case where we do not have preprocessing access to $A$, the vector $x$ is all ones, and $A$ is all $1$s except for one randomly selected entry which has value $2$.
In this case, any algorithm for PAC computation of softmax$(Ax)$ with $\delta < 1-1/n$ (even just identification of the largest entry of $Ax$) requires $\Omega(nd)$ samples.
More practically though, these vectors are the result of a machine learning pipeline, and not of adversarial construction.
As shown by our simulations, this worst case scenario never occurs in practice, and arm pulls are generally well approximated by a Gaussian (see Figure 1 in the attached pdf, now added to the Appendix).
Additionally, note that for any fixed problem instance (fixed matrix $A$ and $x$), all arm pulls are bounded, and are thus sub-Gaussian.
We have added this more detailed discussion to the Appendix, and have added a reference to it in the main text.

2. **Comparison:** While many algorithms have been devised for accelerating softmax computation, no existing methods, to our knowledge, provide $(\epsilon,\delta)$-PAC guarantees (as we formulate in equation (4)) save for exact computation.
Current approximations use some combination of Locality Sensitive Hashing (LSH) to cluster the vocabulary, truncation of the tail to approximate the normalization factor, and sketching to approximate the matrix multiplication.
However, these methods are not truly instance adaptive, and fall prey to many common flaws in machine learning pipelines. For example, reviewer wft8 references a method which learns a quantization of the matrix $A$ during training [2]. The method proposed in this paper has no theoretical guarantees, requires prior knowledge of the correct output label, and only adapts to the training data, not the actual data presented at inference. Since it focuses on learning a good quantization over the training data, it can suffer from potential distributional shifts of the training and test datasets (concerns raised by reviewer kEa9, that our method avoids). See $\textbf{Comparison}$ response to Reviewer kEa9 for more details.


3. **Converting FLOP gains to wall clock gains:** The focus of this paper is to develop the first provably adaptive softmax algorithm with PAC guarantees, highlighting its dramatically improved sample complexity across a wide variety of models and tasks. 
The eventual goal of this method is to enable wall-clock improvements in hardware implementations.
We provided a brief discussion of this in the Limitations and future work section (Lines 333-352), but have now added additional discussion (below) there and to the supplement.
These next steps of converting our provably and empirically performant method into a hardware optimized wall-clock efficient algorithm is an exciting direction of future work, which we detail below.
In most modern-day transformer architectures, memory I/O serves as the primary bottleneck [1].
AdaptiveSoftmax already presents an opportunity to significantly scale down the number of entries of the matrix $A$ that must be loaded at inference time, and, in the future - if memory remains the bottleneck - improve model bandwidth by a similar factor.
This objective appears in reach, since we have designed the components of AdaptiveSoftmax to be amenable to tiling and parallelization.
Most notably, our implementation of AdaptiveSoftmax uses the same column to generate an importance-weighted sample for each active arm.
The reasons for this implementation decision are two-fold.
First, it takes advantage of the locality of entries in the same column to load samples faster, and, second, it removes intra-column correlation, which can yield theoretically improved performance [Baharav and Tse 2019].
Adjacent column samples can also be combined by simply summing their respective importance weights - admitting a simple tiling of our matrix $A$ that could easily be sized particularly to fit individual tiles into SRAM on a GPU along with a copy of the vector $x$ and the current mean/variance estimates for each arm.
Then, we can dynamically load these tiles into SRAM based on the arm estimates as we do currently.
The successive elimination bandits algorithm utilized by AdaptiveSoftmax is also, by choice, quite easily parallelizable.
We may also store two copies of our matrix $A$ — one with wider tiles and one with taller tiles — to take advantage of our tiling at all stages of the AdaptiveSoftmax algorithm: both when a larger number of samples is necessary for fewer arms, in later stages of adaptive estimation, and when a smaller number of samples is necessary for many arms, in earlier stages of the adaptive estimation.
This said, we observe in our experiments that the bulk of compute is invested in our early samples of many arms.
Just using basic parallelization to speed up this step could result in the desired speed improvements.

[1] Ivanov et al. "Data movement is all you need: A case study on optimizing transformers" 2021

[2] Joulin et al. "Efficient softmax approximation for GPUs." 2017.

Meta-review: Computing full softmax can be expensive. This paper proposes a new algorithm to efficiently estimate the top-k softmax values. The algorithm is inspired by the multi-armed bandit problem: it adaptively allocates resources to important output values. The work provides PAC guarantees and empirical results. 

The final recommendations from the reviewers are unanimously acceptance (7, 5, 5, 7), albeit with varying degrees.

The AC agrees with the reviewers that the problem studied is important, and the proposed method is novel and theoretically sound. 

One major concern from the reviewers with less favorable rating is regarding the sub-Gaussianity assumptions. The AC agrees with the authors that this assumption is not strict.

On the other hand, the AC finds that the concerns regarding more discussions and comparisons with other method to be valid (qEfP, wft8).

The topic studied is highly related to sampled softmax. From this point of view, the literature review is not complete . 
The original/early papers discussing sampled softmax are not cited:

*Bengio and Senecal, Quick Training of Probabilistic Neural Nets by Importance Sampling 2003*

*Bengio and Senecal, Adaptive Importance Sampling to Accelerate Training of a Neural Probabilistic Language Model 2008*

There has been many ways of improving the sampling distribution (in the sense of importance sampling framework) since Bengio and Senegal 2008: different types of hashing methods, kernel-based sampling etc. for example the following work and the works it reviewed.

*Blanc and Rendle Adaptive Sampled Softmax with Kernel Based Sampling 2017*

The claim that “no other methods provide theoretical guarantees“ seems to be overstated. Sampled softmax comes with theoretical guarantees from importance sampling directly. See

*Rawat et al Sampled Softmax with Random Fourier Features 2019*

We therefore suggest that the authors work on a better and more comprehensive literature review and slightly tone down the claims in the final version. 

Overall, the AC finds that this is an interesting work to be shared with the community with the above revisions.
Decision: Accept (poster)