Bayesian Coreset Construction via Greedy Iterative Geodesic Ascent

Introduction

Bayesian methods provide a wealth of options for principled parameter estimation and uncertainty quantification. But Markov chain Monte Carlo (MCMC) methods (Robert & Casella, 2004; Neal, 2011; Hoffman & Gelman, 2014), the gold standard for Bayesian inference, typically have complexity $\Theta(NT)$ for dataset size $N$ and number of samples $T$ and are intractable for modern large-scale datasets. Scalable methods (see (Angelino et al., 2016) for a recent survey), on the other hand, often sacrifice the strong guarantees of MCMC and provide unreliable posterior approximations. For example, variational methods and their scalable and streaming variants (Jordan et al., 1999; Wainwright & Jordan, 2008; Hoffman et al., 2013; Ranganath et al., 2014; Broderick et al., 2013; Campbell & How, 2014; Campbell et al., 2015; Dieng et al., 2017) are both susceptible to finding bad local optima in the variational objective and tend to either over- or underestimate posterior variance depending on the chosen discrepancy and variational family.

Bayesian coresets (Huggins et al., 2016; Campbell & Broderick, 2017) provide an alternative approach—based on the observation that large datasets often contain redundant data—in which a small subset of the data of size $M\ll\min\{N,T\}$ is selected and reweighted such that it preserves the statistical properties of the full dataset. The coreset can be passed to a standard MCMC algorithm, providing posterior inference with theoretical guarantees at a significantly reduced $O(M(N+T))$ computational cost. But despite their advantages, existing Bayesian coreset constructions—like many other scalable inference methods—tend to underestimate posterior variance (Fig. 1). This effect is particularly evident when the coreset is small, which is the regime we are interested in for scalable inference.

In this work, we show that existing Bayesian coreset constructions underestimate posterior uncertainty because they scale the coreset log-likelihood suboptimally in order to remain unbiased (Huggins et al., 2016) or to keep their weights in a particular constraint polytope (Campbell & Broderick, 2017). The result is an overweighted coreset with too much “artificial data,” and therefore an overly certain posterior. Taking this intuition to its limit, we demonstrate that there exist models for which previous algorithms output coresets with arbitrarily large relative posterior approximation error at any coreset size (Proposition 2.1). We address this issue by developing a novel coreset construction algorithm, greedy iterative geodesic ascent (GIGA), that optimally scales the coreset log-likelihood to best fit the full dataset log-likelihood. GIGA has the same computational complexity as the current state of the art, but its optimal log-likelihood scaling leads to uniformly bounded relative error for all models, as well as asymptotic exponential error decay (Theorems 3.1 and 3.2). The paper concludes with experimental validation of GIGA on a synthetic vector approximation problem as well as regression models applied to multiple real and synthetic datasets.

Bayesian Coresets

In Bayesian statistical modeling, we are given a dataset $\left(y_{n}\right)_{n=1}^{N}$ of $N$ observations, a likelihood $p(y_{n}|\theta)$ for each observation given a parameter $\theta\in\Theta$, and a prior density $\pi_{0}(\theta)$ on $\Theta$. We assume that the data are conditionally independent given $\theta$. The Bayesian posterior is given by

where the log-likelihood $\mathcal{L}(\theta)$ is defined by

and $Z$ is the marginal likelihood. MCMC returns approximate samples from the posterior, which can be used to construct an empirical approximation to the posterior distribution. Since each sample requires at least one full likelihood evaluation—typically an $\Theta(N)$ operation—MCMC has $\Theta(NT)$ complexity for $T$ posterior samples.

Solving Eq. 3 exactly is not tractable for large $N$ due to the cardinality constraint; approximation is required. Given a norm induced by an inner product, and defining

Campbell & Broderick (2017) replace the cardinality constraint in Eq. 3 with a simplex constraint,

Eq. 5 can be solved while ensuring $\|w\|_{0}\leq M$ using either importance sampling (IS) or Frank–Wolfe (FW) (Frank & Wolfe, 1956). Both procedures add one data point to the linear combination $\mathcal{L}(w)$ at each iteration; IS chooses the new data point i.i.d. with probability $\nicefrac{{\sigma_{n}}}{{\sigma}}$, while FW chooses the point most aligned with the residual error. This difference in how the coreset is built results in different convergence behavior: FW exhibits geometric convergence $\left\|\mathcal{L}(w)-\mathcal{L}\right\|^{2}=O(\nu^{M})$ for some $0<\nu<1$, while IS is limited by the Monte Carlo rate $\left\|\mathcal{L}(w)-\mathcal{L}\right\|^{2}=O(M^{-1})$ with high probability (Campbell & Broderick, 2017, Theorems 4.1, 4.4). For this reason, FW is the preferred method for coreset construction.

Eq. 3 is a special case of the sparse vector approximation problem, which has been studied extensively in past literature. Convex optimization formulations—e.g. basis pursuit (Chen et al., 1999), LASSO (Tibshirani, 1996), the Dantzig selector (Candès & Tao, 2007), and compressed sensing (Candès & Tao, 2005; Donoho, 2006; Boche et al., 2015)—are expensive to solve compared to our greedy approach, and often require tuning regularization coefficients and thresholding to ensure cardinality constraint feasibility. Previous greedy iterative algorithms—e.g. (orthogonal) matching pursuit (Mallat & Zhang, 1993; Chen et al., 1989; Tropp, 2004), Frank–Wolfe and its variants (Frank & Wolfe, 1956; Guélat & Marcotte, 1986; Jaggi, 2013; Lacoste-Julien & Jaggi, 2015; Locatello et al., 2017), Hilbert space vector approximation methods (Barron et al., 2008), kernel herding (Chen et al., 2010), and AdaBoost (Freund & Schapire, 1997)—have sublinear error convergence unless computationally expensive correction steps are included. In contrast, the algorithm developed in this paper has no correction steps, no tuning parameters, and geometric error convergence.

2 Posterior Uncertainty Underestimation

The new $\sum_{n}\sigma_{n}w_{n}=\sigma$ constraint in Eq. 5 has an unfortunate practical consequence: both IS and FW must scale the coreset log-likelihood suboptimally—roughly, by $\sigma$ rather than $\|\mathcal{L}\|$ as they should—in order to maintain feasibility. Since $\sigma\geq\|\mathcal{L}\|$, intuitively the coreset construction algorithms are adding too much “artificial data” via the coreset weights, resulting in an overly certain posterior approximation. It is worth noting that this effect is apparent in the error bounds developed by Campbell & Broderick (2017), which are all proportional to $\sigma$ rather than $\|\mathcal{L}\|$ as one might hope for when approximating $\mathcal{L}$.

which can be made as large as desired by increasing $N$. The result follows since both FW and IS generate a feasible solution $w$ for Eq. 5 satisfying $\|w\|_{0}\leq M$. ∎

In contrast, the red coreset approximation obtained by solving Eq. 3 scales its weight to minimize $\|\mathcal{L}(w)-\mathcal{L}\|$, resulting in a significantly better approximation of posterior uncertainty. Note that we can scale the weight vector $w$ by any $\alpha\geq 0$ without affecting feasibility in the cardinality constraint, i.e., $\|\alpha w\|_{0}\leq\|w\|_{0}$. In the following section, we use this property to develop a greedy coreset construction algorithm that, unlike FW and IS, maintains optimal log-likelihood scaling in each iteration.

Greedy Iterative Geodesic Ascent (GIGA)

In this section, we provide a new algorithm for Bayesian coreset construction and demonstrate that it yields improved approximation error guarantees proportional to $\|\mathcal{L}\|$ (rather than $\sigma\geq\|\mathcal{L}\|$). We begin in Section 3.1 by solving for the optimal log-likelihood scaling analytically. After solving this “radial optimization problem,” we are left with a new optimization problem on the unit hyperspherical manifold. In Sections 3.2 and 3.3, we demonstrate how to solve this new problem by iteratively building the coreset one point at a time. Since this procedure selects the point greedily based on a geodesic alignment criterion, we call it greedy iterative geodesic ascent (GIGA), detailed in Algorithm 1. In Section 3.4 we scale the resulting coreset optimally using the procedure developed in Section 3.1. Finally, in Section 3.5, we prove that Algorithm 1 provides approximation error that is proportional to $\|\mathcal{L}\|$ and geometrically decaying in $M$, as shown by Theorem 3.1.

The particulars of the sequence $\left(\nu_{M}\right)_{M=1}^{\infty}$ are somewhat involved; see Section 3.5 for the detailed development. A straightforward consequence of Theorem 3.1 is that the issue described in Proposition 2.1 has been resolved; the solution $\mathcal{L}(w)$ is always scaled optimally, leading to decaying relative error. Corollary 3.2 makes this notion precise.

For any set of vectors $\left(\mathcal{L}_{n}\right)_{n=1}^{N}$, Algorithm 1 provides a solution to Eq. 3 with error $\leq 1$ relative to $\|\mathcal{L}\|$.

We begin again with the problem of coreset construction for a collection of vectors $\left(\mathcal{L}_{n}\right)_{n=1}^{N}$ given by Eq. 3. Without loss of generality, we assume $\|\mathcal{L}\|>0$ and $\forall n\in[N]$, $\|\mathcal{L}_{n}\|>0$; if $\|\mathcal{L}\|=0$ then $w=0$ is a trivial optimal solution, and any $\mathcal{L}_{n}$ for which $\|\mathcal{L}_{n}\|=0$ can be removed from the coreset without affecting the objective in Eq. 3. Recalling that the weights of a coreset can be scaled by an arbitrary constant $\alpha\geq 0$ without affecting feasibility, we rewrite Eq. 3 as

Taking advantage of the fact that $\|\cdot\|$ is induced by an inner product, we can expand the objective in Eq. 9 as a quadratic function of $\alpha$, and analytically solve the optimization in $\alpha$ to yield

In other words, $\mathcal{L}(w)$ should be rescaled to have norm $\|\mathcal{L}\|$, and then scaled further depending on its directional alignment with $\mathcal{L}$. We substitute this result back into Eq. 9 to find that coreset construction is equivalent to solving

2 Coreset Initialization

For any $\xi\in\Delta^{N-1}$, we have that

By assumption (see Section 3.1), we have $\|\mathcal{L}\|>0$. ∎

3 Greedy Point Selection and Reweighting

and we select the data point at index $n_{t}$ where

Constraining $\gamma\in$ ensures that the resulting $w_{t+1}$ is feasible for Eq. 12. Taking the derivative and setting it to 0 yields the unconstrained optimum of Eq. 19,

It is sufficient to show that both $\zeta_{0}-\zeta_{1}\zeta_{2}$ and $\zeta_{1}-\zeta_{0}\zeta_{2}$ are nonnegative and that at least one is strictly positive. First, examining Eq. 18, note that for any $\xi\in\Delta^{N-1}$,

4 Output

5 Convergence Guarantee

The constants $\tau$ and $\epsilon$ satisfy

The geodesic alignment $\left<d_{t},d_{tn_{t}}\right>$ satisfies

The proof of Theorem 3.1 below follows by using the \tau\mathchoice{{\hbox{\displaystyle\sqrt{J_{t}\,}}\lower 0.4pt\hbox{\vrule height=6.83331pt,depth=-5.46667pt}}}{{\hbox{\textstyle\sqrt{J_{t}\,}}\lower 0.4pt\hbox{\vrule height=6.83331pt,depth=-5.46667pt}}}{{\hbox{\scriptstyle\sqrt{J_{t}\,}}\lower 0.4pt\hbox{\vrule height=4.78333pt,depth=-3.82668pt}}}{{\hbox{\scriptscriptstyle\sqrt{J_{t}\,}}\lower 0.4pt\hbox{\vrule height=3.41666pt,depth=-2.73334pt}}} bound from Lemma 3.6 to obtain an $O(1/t)$ bound on $J_{t}$, and then combining that result with the $f(J_{t})$ bound from Lemma 3.6 to obtain the desired geometric decay.

Substituting the \tau\mathchoice{{\hbox{\displaystyle\sqrt{J_{t}\,}}\lower 0.4pt\hbox{\vrule height=6.83331pt,depth=-5.46667pt}}}{{\hbox{\textstyle\sqrt{J_{t}\,}}\lower 0.4pt\hbox{\vrule height=6.83331pt,depth=-5.46667pt}}}{{\hbox{\scriptstyle\sqrt{J_{t}\,}}\lower 0.4pt\hbox{\vrule height=4.78333pt,depth=-3.82668pt}}}{{\hbox{\scriptscriptstyle\sqrt{J_{t}\,}}\lower 0.4pt\hbox{\vrule height=3.41666pt,depth=-2.73334pt}}} bound from Lemma 3.6 into Eq. 28 and applying a standard inductive argument (e.g. the proof of Campbell & Broderick (2017, Lemma A.6)) yields

Multiplying by $\|\mathcal{L}\|^{2}$, taking the square root, and noting that $\eta$ from Eq. 8 is equal to $\textstyle\sqrt{J_{1}\,}$ gives the final result. The convergence of $f(B(t))\to\epsilon$ as $t\to\infty$ shows that the rate of decay in the theorem statement is \nu=\mathchoice{{\hbox{\displaystyle\sqrt{1-\epsilon^{2}\,}}\lower 0.4pt\hbox{\vrule height=6.44444pt,depth=-5.15558pt}}}{{\hbox{\textstyle\sqrt{1-\epsilon^{2}\,}}\lower 0.4pt\hbox{\vrule height=6.44444pt,depth=-5.15558pt}}}{{\hbox{\scriptstyle\sqrt{1-\epsilon^{2}\,}}\lower 0.4pt\hbox{\vrule height=4.51111pt,depth=-3.6089pt}}}{{\hbox{\scriptscriptstyle\sqrt{1-\epsilon^{2}\,}}\lower 0.4pt\hbox{\vrule height=3.44165pt,depth=-2.75334pt}}}. ∎

Experiments

In this section we evaluate the performance of GIGA coreset construction compared with both uniformly random subsampling and the Frank–Wolfe-based method of Campbell & Broderick (2017). We first test the algorithms on simple synthetic examples, and then test them on logistic and Poisson regression models applied to numerous real and synthetic datasets. Code for these experiments is available at https://github.com/trevorcampbell/bayesian-coresets.

To generate Fig. 1, we generated $\mu\sim\mathcal{N}(0,1)$ followed by $(y_{n})_{n=1}^{10}\sim\mathcal{N}(\mu,1)$. We then constructed Bayesian coresets via FW and GIGA using the norm specified in (Campbell & Broderick, 2017, Section 6). Across 1000 replications of this experiment, the median relative error in posterior variance approximation is 3% for GIGA and 48% for FW.

2 Synthetic Vector Sum Approximation

In this experiment, we generated 20 independent datasets consisting of $10^{6}$ 50-dimensional vectors from the multivariate normal distribution with mean 0 and identity covariance. We then constructed coresets for each of the datasets via uniformly random subsampling (RND), Frank–Wolfe (FW), and GIGA. We compared the algorithms on two metrics: reconstruction error, as measured by the 2-norm between $\mathcal{L}(w)$ and $\mathcal{L}$; and representation efficiency, as measured by the size of the coreset.

Fig. 3 shows the results of the experiment, with reconstruction error in Fig. 3(a) and coreset size in Fig. 3(b). Across all construction iterations, GIGA provides a 2–4 order-of-magnitude reduction in error as compared with FW, and significantly outperforms RND. The exponential convergence of GIGA is evident. The flat section of FW/GIGA in Fig. 3(a) for iterations beyond $10^{2}$ is due to the algorithm reaching the limits of numerical precision. In addition, Fig. 3(b) shows that GIGA can improve representational efficiency over FW, ceasing to grow the coreset once it reaches a size of 120, while FW continues to add points until it is over twice as large. Note that this experiment is designed to highlight the strengths of GIGA: by the law of large numbers, as $N\to\infty$ the sum $\mathcal{L}$ of the $N$ i.i.d. standard multivariate normal data vectors satisfies $\|\mathcal{L}\|\to 0$, while the sum of their norms $\sigma\to\infty$ a.s. But Section A.1 shows that even in pathological cases, GIGA outperforms FW due to its optimal log-likelihood scaling.

3 Bayesian Posterior Approximation

We constructed coresets for each of the datasets via uniformly random subsampling (RND), Frank–Wolfe (FW), and GIGA using the weighted Fisher information distance,

where $\mathcal{L}_{n}(\theta)$ is the log-likelihood for data point $n$, and $\hat{\pi}$ is obtained via the Laplace approximation. In order to do so, we approximated all vectors $\mathcal{L}_{n}$ using a 500-dimensional random feature projection (following Campbell & Broderick, 2017). For posterior inference, we used Hamiltonian Monte Carlo (Neal, 2011) with 15 leapfrog steps per sample. We simulated a total of 6,000 steps, with 1,000 warmup steps for step size adaptation with a target acceptance rate of 0.8, and 5,000 posterior sampling steps. Appendix A shows results for similar experiments using random-walk Metropolis–Hastings and the No-U-Turn Sampler (Hoffman & Gelman, 2014). We ran 20 trials of projection / coreset construction / MCMC for each combination of dataset, model, and algorithm. We evaluated the coresets at logarithmically-spaced construction iterations between $M=1$ and $1000$ by their median posterior Fisher information distance, estimated using samples obtained by running posterior inference on the full dataset in each trial. We compared these results versus computation time as well as coreset size. The latter serves as an implementation-independent measure of total cost, since the cost of running MCMC is much greater than running coreset construction and depends linearly on coreset size.

The results of this experiment are shown in Fig. 4. In Fig. 4(a) and Fig. 4(b), the Fisher information distance is normalized by the median distance of RND for comparison. In Fig. 4(b), the computation time is normalized by the median time to run MCMC on the full dataset. The suboptimal coreset log-likelihood scaling of FW can be seen in Fig. 4(a) for small coresets, resulting in similar performance to RND. In contrast, GIGA correctly scales posterior uncertainty across all coreset sizes, resulting in a major (3–4 orders of magnitude) reduction in error. Fig. 4(b) shows the same results plotted versus total computation time. This confirms that across a variety of models and datasets, GIGA provides significant improvements in posterior error over the state of the art.

Conclusion

This paper presented greedy iterative geodesic ascent (GIGA), a novel Bayesian coreset construction algorithm. Like previous algorithms, GIGA is simple to implement, has low computational cost, and has no tuning parameters. But in contrast to past work, GIGA scales the coreset log-likelihood optimally, providing significant improvements in the quality of posterior approximation. The theoretical guarantees and experimental results presented in this work reflect this improvement.

Acknowledgments

This research is supported by an MIT Lincoln Laboratory Advanced Concepts Committee Award, ONR grant N00014-17-1-2072, a Google Faculty Research Award, and an ARO YIP Award.

References

Appendix A Additional results

A.2 Alternate inference algorithms

We reran the same experiment as described in Section 4.3, except we swapped the inference algorithm for random-walk Metropolis–Hastings (RWMH) and the No-U-Turn Sampler (NUTS) (Hoffman & Gelman, 2014). When using RWMH, we simulated a total of 50,000 steps: 25,000 warmup steps including covariance adaptation with a target acceptance rate of 0.234, and 25,000 sampling steps thinned by a factor of 5, yielding 5,000 posterior samples. If the acceptance rate for the latter 25,000 steps was not between 0.15 and 0.7, we reran the procedure. When using NUTS, we simulated a total of 6,000 steps: 1,000 warmup steps including leapfrog step size adaptation with a target acceptance rate of 0.8, and 5,000 sampling steps.

The results for these experiments are shown in Figs. 6 and 7, and generally corroborate the results from the experiments using Hamiltonian Monte Carlo in the main text. One difference when using NUTS is that the performance versus computation time appears to follow an “S”-shaped curve, which is caused by the dynamic path-length adaptation provided by NUTS. Consider the log-likelihood of logistic regression, which has a “nearly linear” region and a “nearly flat” region. When the coreset is small, there are directions in latent space that point along “nearly flat” regions; along these directions, u-turns happen only after long periods of travel. When the coreset reaches a certain size, these “nearly flat” directions are all removed, and u-turns happen more frequently. Thus we expect the computation time as a function of coreset size to initially increase smoothly, then drop quickly, followed by a final smooth increase, in agreement with Fig. 7(b).

Appendix B Technical Results and Proofs

We begin with the \tau\mathchoice{{\hbox{\displaystyle\sqrt{J_{t}\,}}\lower 0.4pt\hbox{\vrule height=6.83331pt,depth=-5.46667pt}}}{{\hbox{\textstyle\sqrt{J_{t}\,}}\lower 0.4pt\hbox{\vrule height=6.83331pt,depth=-5.46667pt}}}{{\hbox{\scriptstyle\sqrt{J_{t}\,}}\lower 0.4pt\hbox{\vrule height=4.78333pt,depth=-3.82668pt}}}{{\hbox{\scriptscriptstyle\sqrt{J_{t}\,}}\lower 0.4pt\hbox{\vrule height=3.41666pt,depth=-2.73334pt}}} bound. For any $\xi\in\Delta^{N-1}$,

Maximizing over all valid choices of $s$ yields

Next, we develop the $f(J_{t})$ bound. Note that

We add $\{0\}$ into the minimization since $\beta\leq 0$ guarantees that the derivative of the above with respect to $J_{t}$ is nonpositive (which we will require in proving the main theorem). For all $J_{t}$ small enough such that \mathchoice{{\hbox{\displaystyle\sqrt{1-J_{t}\,}}\lower 0.4pt\hbox{\vrule height=6.83331pt,depth=-5.46667pt}}}{{\hbox{\textstyle\sqrt{1-J_{t}\,}}\lower 0.4pt\hbox{\vrule height=6.83331pt,depth=-5.46667pt}}}{{\hbox{\scriptstyle\sqrt{1-J_{t}\,}}\lower 0.4pt\hbox{\vrule height=4.78333pt,depth=-3.82668pt}}}{{\hbox{\scriptscriptstyle\sqrt{1-J_{t}\,}}\lower 0.4pt\hbox{\vrule height=3.41666pt,depth=-2.73334pt}}}\mathchoice{{\hbox{\displaystyle\sqrt{1-\beta^{2}\,}}\lower 0.4pt\hbox{\vrule height=8.74889pt,depth=-6.99915pt}}}{{\hbox{\textstyle\sqrt{1-\beta^{2}\,}}\lower 0.4pt\hbox{\vrule height=8.74889pt,depth=-6.99915pt}}}{{\hbox{\scriptstyle\sqrt{1-\beta^{2}\,}}\lower 0.4pt\hbox{\vrule height=6.14998pt,depth=-4.92001pt}}}{{\hbox{\scriptscriptstyle\sqrt{1-\beta^{2}\,}}\lower 0.4pt\hbox{\vrule height=4.7611pt,depth=-3.8089pt}}}\epsilon+\mathchoice{{\hbox{\displaystyle\sqrt{J_{t}\,}}\lower 0.4pt\hbox{\vrule height=6.83331pt,depth=-5.46667pt}}}{{\hbox{\textstyle\sqrt{J_{t}\,}}\lower 0.4pt\hbox{\vrule height=6.83331pt,depth=-5.46667pt}}}{{\hbox{\scriptstyle\sqrt{J_{t}\,}}\lower 0.4pt\hbox{\vrule height=4.78333pt,depth=-3.82668pt}}}{{\hbox{\scriptscriptstyle\sqrt{J_{t}\,}}\lower 0.4pt\hbox{\vrule height=3.41666pt,depth=-2.73334pt}}}\beta\geq 0, the derivative of the above with respect to $\left<-d,d_{\infty n}\right>$ is nonnegative. Therefore, minimizing yields

which holds for any such small enough $J_{t}$. But note that we’ve already proven the \left<d_{t},d_{tn_{t}}\right>\geq\tau\mathchoice{{\hbox{\displaystyle\sqrt{J_{t}\,}}\lower 0.4pt\hbox{\vrule height=6.83331pt,depth=-5.46667pt}}}{{\hbox{\textstyle\sqrt{J_{t}\,}}\lower 0.4pt\hbox{\vrule height=6.83331pt,depth=-5.46667pt}}}{{\hbox{\scriptstyle\sqrt{J_{t}\,}}\lower 0.4pt\hbox{\vrule height=4.78333pt,depth=-3.82668pt}}}{{\hbox{\scriptscriptstyle\sqrt{J_{t}\,}}\lower 0.4pt\hbox{\vrule height=3.41666pt,depth=-2.73334pt}}} bound, which is always nonnegative; so the only time the current bound is “active” is when it is itself nonnegative, i.e. when $J_{t}$ is small enough. Therefore the bound

Appendix C Cap-tree Search

When choosing the next point to add to the coreset, we need to solve the following maximization with $O(N)$ complexity:

If we write $\zeta=\alpha_{u}u+\alpha_{v}v+\sum_{i}\alpha_{i}z_{i}$ where $z_{i}$ completes the basis of $u,v$ etc, and $\xi=\beta_{u}u+\beta_{v}v+\sum_{i}\beta_{i}z_{i}$,

Noting that $\alpha_{i}$ doesn’t appear in the objective, we maximize $\sum_{i}\alpha_{i}\beta_{i}$ to find the equivalent optimization

where the norm on $|\beta_{v}|$ comes from the fact that we can choose the sign of $\alpha_{v}$ arbitrarily, ensuring the optimum has $\alpha_{v}\geq 0$. Now define

Since $\eta$ is now decoupled from the optimization, we can solve

to make the feasible region in $\gamma$ as large as possible. If $|\beta_{v}|>r$, we maximize Eq. 71 by sending $\eta\to\infty$ yielding a maximum of 1 in the original optimization. Otherwise, note that at $\eta=1$ the derivative of the objective is $+\infty$, so we know the constraint $\eta=1$ is not active. Therefore, taking the derivative and setting it to 0 yields

Substituting back into the original optimization,

If \beta_{u}\geq\mathchoice{{\hbox{\displaystyle\sqrt{r^{2}-|\beta_{v}|^{2}\,}}\lower 0.4pt\hbox{\vrule height=9.30444pt,depth=-7.44359pt}}}{{\hbox{\textstyle\sqrt{r^{2}-|\beta_{v}|^{2}\,}}\lower 0.4pt\hbox{\vrule height=9.30444pt,depth=-7.44359pt}}}{{\hbox{\scriptstyle\sqrt{r^{2}-|\beta_{v}|^{2}\,}}\lower 0.4pt\hbox{\vrule height=6.53888pt,depth=-5.23112pt}}}{{\hbox{\scriptscriptstyle\sqrt{r^{2}-|\beta_{v}|^{2}\,}}\lower 0.4pt\hbox{\vrule height=5.03888pt,depth=-4.03113pt}}}, then $\gamma=1$ is feasible and the optimum is 1. Otherwise, note that at $\gamma=-1$, the derivative of the constraint is $+\infty$ and the derivative of the objective is $1$, so the constraint $\gamma=-1$ is not active. Therefore, we can solve the unconstrained optimization by taking the derivative and setting to 0, yielding

Therefore, the upper bound is as follows:

Appendix D Datasets

The Phishing dataset is available online at https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/binary.html. The DS1 dataset is available online at http://komarix.org/ac/ds/. The BikeTrips dataset is available online at http://archive.ics.uci.edu/ml/datasets/Bike+Sharing+Dataset. The AirportDelays dataset was constructed using flight delay data from http://stat-computing.org/dataexpo/2009/the-data.html and historical weather information from https://www.wunderground.com/history/.