Mixability made efficient: Fast online multiclass logistic regression
Year: 2022

Abstract:
Mixability has been shown to be a powerful tool to obtain algorithms with optimal regret. However, the resulting methods often suffer from high computational complexity which has reduced their practical applicability. For example, in the case of multiclass logistic regression, the aggregating forecaster (Foster et al. (2018)) achieves a regret of O(log(Bn)) whereas Online Newton Step achieves O(e log(n)) obtaining a double exponential gain in B (a bound on the norm of comparative functions). However, this high statistical performance is at the price of a prohibitive computational complexity O(n). In this paper, we use quadratic surrogates to make aggregating forecasters more efficient. We show that the resulting algorithm has still high statistical performance for a large class of losses. In particular, we derive an algorithm for multi-class logistic regression with a regret bounded by O(B log(n)) and a computational complexity of only O(n).

heading: 1 Introduction
text: In online learning, a learner sequentially interacts with an environment and tries to learn based on data observed on the fly (Cesa-Bianchi and Lugosi (2006), Hazan et al. (2016)). More formally, at each iteration t ≥ 1, the learner receives an input xt in some space X ; makes a prediction ŷt in a decision domain Ŷ and the environment reveals the output yt ∈ Y . The inputs xt and the outputs yt are sequentially chosen by the environment and can be arbitrary. No stochastic assumption (except boundedness) on the data sequence (xt, yt)1≤t≤n is made. The accuracy of a prediction ŷt ∈ Ŷ at instant t ≥ 1 for the outcome yt ∈ Y is measured through a loss function ` : Ŷ × Y → R. The learner aims at minimizing his cumulative regret
Rn(f) = n∑ t=1 ` ( ŷt, yt ) − n∑ t=1 ` ( f(xt), yt ) , (1)
with respect to any function f ∈ F in a reference class. All along this paper, we will consider parametric class of functions fθ indexed by θ ∈ Θ ⊂ Rd. We will also assume convexity of the loss according to index θ, which we denote `x,y : θ 7→ `(fθ(x), y) for any (x, y) ∈ X × Y . In this
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
general context, many algorithms have been designed based on different assumptions and obtaining different trade-offs for regret. We review below the most relevant ones for our purpose.
Assuming only convexity of the loss `, a well-known strategy for the learner is Online Gradient Descent (Zinkevich (2003)) with an optimal regret of orderO( √ n). However, if the loss is α-mixable, the learner may achieve the faster rate O( 1α log(n)) by using an aggregating forecaster (see Vovk (2001), Van Erven et al. (2015)). Yet, such an algorithm is not constructive in general and when it is, the computational complexity is often very high (a notable exception is the least-squares setting). To reduce it, a stronger assumption has been introduced by Hazan et al. (2007) with η-exp concavity. Under the latter hypothesis, Online Newton Step (ONS) has a regret bounded by O( 1η log(n)) and a computational complexity of O(n). However this efficiency comes at the price of deteriorating statistical performance. Indeed, η-exp concavity implies α-mixability for η ≥ α and in some cases, the gap between η and α can be very large.
The most spectacular case of this phenomenon occurs for logistic regression. In this setting, the loss is defined as
` (ŷt, yt) = − log (σ(ŷt)yt) where σ(z)i = ezi∑ j e zj
and the regret is computed with respect to linear functions F = {x 7→ Wx,W ∈ Θ}, where Θ = B(RK×d, B) is the Frobenius bounded ball of radius B > 0. On this subset, the logistic loss is η-exp concave only for η ≤ e−B . The regret of Online Newton Step can thus be of order O(eB log(n)). On the other hand, as remarked by Kakade and Ng (2008) (binary case) and Foster et al. (2018), the logistic loss is 1-mixable and an aggregating algorithm may achieve O(log(nB)). Nevertheless, the algorithm relies on Monte Carlo methods and must sample from log-concave probability distributions which is extremely computationally expensive (see Table 1). Therefore, in this framework, the exp-concavity assumption leads to an efficient algorithm with low dependence on B while mixability yields an inefficient algorithm with much better statistical performance.
Between these two extremes, other algorithms and trade-offs have been analyzed. First, it has been shown that in some situations, Follow The Regularized Leader (FTRL) can achieve a fast rate without an exponential constant (Bach et al. (2010), Marteau-Ferey et al. (2019), Ostrovskii et al. (2021)). However, several additional assumptions are necessary to achieve these rates. Unfortunately, they are essentially unavoidable. Indeed, Hazan et al. (2014) showed a polynomial lower bound for the proper algorithm (i.e., with linear prediction function) in the regime B = log(n). This prevents algorithms like ONS and Follow The Regularized Leader from reaching logarithmic regret without an exponential constant in B. This result motivates the search for improper algorithms with better computational complexity than Foster et al. (2018).
For binary logistic regression, two efficient improper algorithms have been proposed in the literature. First, in the i.i.d. framework, Mourtada and Gaïffas (2019) have proposed the Sample Minmax Predictor (SMP). It achieves an excess risk of order O((B2 + d)/n) with computational complexity equivalent to Follow The Regularized leader. In the online framework, Jézéquel et al. (2020) proposed AIOLI, which is based on quadratic approximations of the logistic loss as well as virtual labels to regularize. The regret is upper-bounded O(dB log(n)) and the computational complexity is O(n(d2 + log(n))).
These previous works left open the question of achieving the same type of performance in a setting other than binary logistic regression. In particular, for multi-class logistic regression, no other algorithm than Foster et al. (2018) is known to achieve logarithmic regret without an exponential constant.
Contributions We introduce in Section 3 a new generic online learning algorithm, that we call Gaussian Aggregating Forecaster (GAF) . GAF achieves logarithmic regret with a small multiplicative constant for a large class of convex loss functions (see Theorem 1). The latter includes several popular loss functions such as squared loss, binary and multi-class logistic loss. Our assumptions on the loss functions are slightly stronger than α-mixability but generally weaker than assumptions widely used in the statistical framework such as generalized self-concordance (see Bach et al. (2010)).
In the particular but significant setting of multi-class logistic regression, GAF has a regret bounded by O(dKB2 + dKB log(n)) and a total computational complexity of O(nK3d2 + K2n4), thus significantly improving on the O(B6n25(Bn + dK)12) complexity of Foster et al. (2018), which was the best known to date for algorithms without exponential dependence on B. Table 1 summarizes the regrets and computational cost obtained by the relevant algorithms in this framework. It is worth pointing out that, by using standard online-to-batch conversion Helmbold and Warmuth (1995), this paper also provides new results on the excess risk in the statistical i.i.d. Indeed, even in this extensively used framework, the only existing improper algorithms without exponential constant in B so far are GAF and Foster et al. (2018). Jézéquel et al. (2020) and Mourtada and Gaïffas (2019) are in fact restricted to binary outputs only.
Our new algorithm is inspired by the aggregating forecaster of Vovk (2001) that we apply to quadratic approximations of the loss, in order to make it more efficient. The high-level idea is that these quadratic approximations replace complicated log-concave distributions with Gaussian distributions, from which it is much easier to sample and thereby drastically reducing the complexity of the algorithm. We believe that this approach can be applied to much broader contexts than those analyzed here.

heading: 2 Setting
text: We recall here the setting and introduce the main notations and assumptions that will be used throughout the paper.
Setting and notation Our framework is formalized as a sequential game between a learner and an environment. At each forecasting instance t ≥ 1, the learner is given an input xt ∈ X ; forms a prediction ŷt ∈ Ŷ ⊆ RK (possibly based on the current input xt and on the past information x1, y1, . . . , xt−1, yt−1). Note that the prediction space Ŷ ⊆ RK may be uni-dimensional (K = 1) in some settings (least-square regression) or multi-dimensional (K ≥ 1) in some cases (e.g., vectorvalued regression, multi-class classification). Then, the environment chooses yt ∈ Y; reveals it to the learner which incurs the loss `(ŷt, yt). Finally, the performance of the learner is assessed by the regret
Rn(θ) = n∑ t=1 `(ŷt, yt)− n∑ t=1 `xt,yt(θ) ,
with respect to all θ ∈ Θ ⊂ Rd. Here, `x,y(θ) = `(fθ(x), y) where F = {fθ : X → Ŷ} is the reference class of functions. For simplicity of notation, we also define `t(θ) = `xt,yt(θ) and Lt(θ) = ∑t s=1 `s(θ) + λ‖θ‖2 for all t ≥ 1 and λ > 0.
In all specific examples considered in this work, the learner will be compared to linear functions fθ : x 7→ θ>Φ(x), where θ ∈ Θ = B(Rd, B) and Φ is a function from X to Rd×K such that for any x ∈ X and i ∈ [K], ‖Φ(x).,i‖ ≤ R. We introduce and discuss below our main assumptions on the losses.
Assumptions We assume that, for all (x, y) ∈ X × Y , the losses `x,y are convex, C2 and satisfy the following assumptions.
(A1) The loss function ` is α-mixable. In other words, for all (Gaussian) probability distributions π over Rd and input x ∈ X , there exists ŷ ∈ Ŷ such that for all y ∈ Y ,
`(ŷ, y) ≤ − 1 α
log ( Eθ∼πe−α`x,y(θ) ) .
(A2) There exists ζ > 0 such that, for all (x, y) ∈ X × Y and θ1, θ2 ∈ Rd,
`x,y(θ1) ≤ `x,y(θ2) +∇`x,y(θ2)>(θ1 − θ2) + eζ‖θ1−θ2‖ 2 ‖θ1 − θ2‖2∇2`x,y(θ2) .
(A3) There exists β > 0 such that, for all (x, y) ∈ X × Y and (θ1, θ2) ∈ Θ× Rd,
`x,y(θ1) ≥ `x,y(θ2) +∇`x,y(θ2)>(θ1 − θ2) + β
2 (θ1 − θ2)>∇2`x,y(θ2)(θ1 − θ2) .
(A4) ` is γ-smooth i.e., for all x ∈ X , y ∈ Y, θ ∈ Rd,∇2`x,y(θ) ≤ γI . In general, these assumptions are rather weak and related to standard assumptions of (online) convex optimization. (A1) follows from instance from exp-concavity. As we illustrate in Section 4, all are satisfied for broadly used loss functions such as squared loss, binary and multi-class logistic loss. We provide more details on these assumptions below.
Mixability (A1) is a standard assumption in online learning to achieve fast rates for the regret (Van Erven et al., 2015), which was introduced by Vovk (2001). It is a weaker condition than exp-concavity because (A1) holds with ŷ = Eθ∼π[θ] for η-exp concave loss functions when η ≥ α. Assumption (A2) basically prevents the Hessian from changing too quickly on θ. This condition is new but not very restrictive because, as we prove in Lemma 6, it is weaker than other widely accepted assumptions such as generalized self-concordance (Bach et al., 2010).
The lower bound (A3) is the most original assumption and is crucial to make aggregating forecaster efficient. This bound is close to the one needed in the analyses of ONS, which can be derived from η-exp-concavity. However, the key difference is that (A3) uses the Hessian in the quadratic term and not the outer product of the gradient. This makes the right-hand-side closer to the third Taylor series approximation of the loss and allows much larger values for β. For example, as we will show in the section 4.2, this change allows us to remove the exponential constant for the logistic regression setting. In this case, (A3) is verified with β ' B−1, while η-exp-concavity would require η ' e−B .

heading: 3 Algorithm and regret bound
text: GAF is a new sequential forecasting rule inspired by the aggregating forecaster of Vovk (2001). It may be implemented if the loss function ` is mixable (A1) and ensures fast performance guarantees on the regret under Assumptions (A1)-(A4) as shown in Theorem 1. GAF requires the following hyper-parameters: a regularization parameter λ > 0 and α, β > 0 such that the loss ` satisfies assumptions (A1)-(A4).
Algorithm (GAF) At each time step, the prediction ŷt is formed by GAF by following a two steps procedure. A first estimator θt ∈ Rd is computed by solving
θt = argmin θ∈Rd { t−2∑ s=1 ˜̀ s(θ) + `t−1(θ) + λ‖θ‖22 } . (2)
Basically, θt follows a regularized leader where the losses `s for s ≤ t − 2 are substituted with quadratic approximations ˜̀s. Then, GAF computes the quadratic approximation of `t−1 at point θt. For any θ ∈ Rd, it defines
˜̀ t−1(θ) = `t−1(θt) +∇`t−1(θt)>(θ − θt) +
β 2 (θ − θt)>∇2`t−1(θt)(θ − θt). (3)
Finally, GAF predicts ŷt ∈ Ŷ such that for all y ∈ Y ,
`(ŷt, y) ≤ − 1 α log ( Eθ∼P̃t−1e −α`xt,y(θ) ) , where P̃t−1(θ) = e−αL̃t−1(θ)∫ Rd e −αL̃t−1(θ)dθ , (4)
with L̃t−1(θ) = ∑t−1 s=1 ˜̀ s(θ) + λ‖θ‖2. Such a prediction ŷt exists as soon as Assumption (A1) is true. In Section 4, we present some specific cases where ŷt can be computed in a closed form.
We now state our main theoretical result, which is an upper bound on the regret suffered by GAF.
Theorem 1. Let d, n ≥ 1, B > 0, and Θ ⊂ B(Rd, B). Let (x1, y1), ..., (xn, yn) ∈ X × Y be an arbitrary sequence of observations and ` a loss function that verifies Assumptions (A1)-(A4) with α, ζ, β, γ > 0. GAF (4), run with regularization parameter λ ≥ max{4, d}ζα−1, satisfies the following upper-bound on the regret
Rn(θ) ≤ λ‖θ‖2 + d
α
[ 1
2 +
2 √ 3
β
] log ( 1 + nβγ
2λ
) , ∀θ ∈ Θ .
In particular, the choice λ = max{4, d}ζα−1 yields
Rn(θ) ≤ max{4, d} ζB2 α + d α
[ 1
2 +
2 √ 3
β
] log ( 1 + nαβγ
dζ
) , ∀θ ∈ Θ .
Theorem 1 states that GAF has a logarithmic regret in the number of samples with a multiplicative constant proportional to α−1(1 + β−1). Recalling that α is the mixability parameter, α−1 is the optimal multiplicative constant that can be obtained using an (computationally expensive) aggregating forecaster. The constant β is the curvature parameter in the quadratic lower-bound assumption (A3). Having β−1 as a multiplicative constant can be seen as similar to the regret bound of ONS. However, it is important to note the crucial difference that, unlike ONS, it is the Hessian that is used in the quadratic substitutes of Assumption (A3) and not the gradient outer product. Thus, in some cases, it is possible to have a much larger value for β, than would be attainable using the exp-concavity. For example, for multiclass logistic regression, the loss is 1-mixable (Proposition 1, Foster et al., 2018) but only e−B-exp concave. We prove in Lemma 4 that this loss verifies Assumption (A3) with parameter β ' (log(K) + BR)−1. Therefore, GAF achieves, in the logistic case, a regret upper bounded by O(dB2R2 + d(log(K) +BR) log(n)). More details on this specific case are given in Section 4.2.
We provide below a sketch of the proof of Theorem 1 that highlights the key steps in the proof. The full proof is available in Appendix B.
Sketch of proof. The proof starts from the definition (4) of the prediction ŷt formed by the algorithm. It satisfies the mixability assumption (A1) for all y ∈ Y
`(ŷt, y) ≤ − 1 α log ( Eθ∼P̃t−1(θ)e −α`t(θ) ) , where P̃t−1(θ) = e−αL̃t−1(θ)∫ Rd e −αL̃t−1(θ)dθ .
Next, we decompose the term on the right side of the inequality into two terms: an approximation error that will be small for well-chosen parameters; and a term, which telescopes and is standard for aggregating forecaster’s analysis (Vovk, 2001), except that here it is written with quadratic surrogate losses. We have,
`(ŷt, yt) ≤ − 1
α log
(∫ Rd e
−α(`t(θ)−˜̀t(θ))−αL̃t(θ)dθ∫ Rd e −αL̃t−1(θ)dθ
)
= − 1 α log ( Eθ∼P̃te −α[`t(θ)−˜̀t(θ)] )
︸ ︷︷ ︸ Ωt
+ 1
α log
(∫ Rd e
−αL̃t−1(θ)dθ∫ Rd e −αL̃t(θ)dθ ) ︸ ︷︷ ︸
Ψt
.
Now, the core idea is that since L̃t is quadratic, the integrals inside Ψt are Gaussian, which allows closed form formulas for the analysis (and fast computational time). Using that ∇˜̀t(θt+1) = ∇`t(θt+1) and the expression of Gaussian integrals, it is possible to show that
Ψt = L̃t(θt+1)− L̃t−1(θt) + 1
2α log ( |At| |At−1| ) ,
where At is the Hessian of L̃t/2. Therefore, summing over t = 1, . . . , n and using L̃0(θ1) = 0, it yields
n∑ t=1 `(ŷt, yt)− L̃n(θn+1) ≤ n∑ t=1 Ωt + 1 2α log ( |An| |A0| ) .
But, from Assumption (A3) together with the definition (2) of θt+1, for any θ ∈ Θ, Lt(θ) ≥ L̃t(θ) ≥ L̃t(θt+1). Thus, the regret can be bounded for any θ ∈ Θ by
Rn(θ) = n∑ t=1 `(ŷt, yt)− Ln(θ) ≤ λ‖θ‖2 + n∑ t=1 Ωt + 1 2α log ( |An| |λI| ) . (5)
Finally, Assumption (A2) allows to bound the approximation error Ωt by a telescopic term close to the usual one but with 1αβ as multiplicative constant instead of 1 α ,
Ωt ≤ C
αβ log ( |At| |At−1| ) .
The proof is concluded by applying Assumption (A4).

heading: 4 Specific settings: squared loss and multi-class logistic
text: In this section we show how our general framework cover several interesting settings. In particular, we prove that our Assumption (A1)-(A4) are satisfied and provide concrete implementations of GAF in those contexts.

heading: 4.1 Squared loss
text: For linear regression with squared loss, we recover classical results by Vovk (2001). This framework is defined, as a special case of our generic setting of Section 2, by setting: input domain X = B(Rd, R), output domainY = [−Y, Y ] with Y > 0, decision domain Ŷ = R, loss function `(ŷt, yt) = (ŷt−yt)2 and input feature Φ(x) = x. Now we can verify that all assumptions are true:
(A1) By Lemma 2 and Lemma 3 of Vovk (2001), the squared loss is mixable with parameter α = 1/(2Y 2). (A2) The inequality is true for any ζ > 0 as ex > 1/2 for all x > 0. We will take ζ = (8Y 2dB2)−1. (A3) It holds with β = 1 because the Taylor expansion of order 2 of a quadratic function is an equality.
(A4) It is true with γ = R2 since ∇2`x,y(θ) = 2xx> ≤ 2R2I . In particular, in this particular case the quadratic approximations are exact ˜̀t = `t for all t ≥ 1. Therefore, GAF reduces to the classical non-linear Ridge regression of Vovk (2001) and Azoury and Warmuth (2001) that predicts
ŷt = θ̂ > t xt with θ̂t = argminθ∈Rd { Lt−1(θ) + θ >xt + λ‖θ‖2 } .
A direct consequence of Theorem 1, by substituting the constants derived above, yields Rn(θ) ≤ 1 + 8Y 2d log ( 1 + 16nB2R2 ) .
This result essentially recovers the existing regret bound for the non-linear Ridge regression (see e.g., Theorem 11.8 of Cesa-Bianchi and Lugosi (2006)).

heading: 4.2 Logistic regression
text: The most interesting application of our general framework is logistic regression. In this section, we show the validity of our assumptions and the concrete implementation of GAF in this context. We present the theoretical results obtained and discuss the computational complexity.
The goal of logistic regression is to form a K-dimensional prediction ŷ ∈ Ŷ = RK of a categorical label y ∈ Y = {1, ...,K} from the observation of an input label x ∈ X = B(Rd′ , R), d′ ≥ 1. The performance of ŷ is measured by the logistic loss defined by
` (ŷ, y) = − log (σ(ŷ)y) where σ(z)i = e zi∑ j e zj .
Defining the input feature Φ(x) ∈ Rd×K with d = d′K and the linear predictions respectively as
Φ(x)i,j = { xi−d′j if d′j ≤ i < d′(j + 1) 0 otherwise
and fθ(x) = θ>Φ(x), for θ ∈ B(Rd, B),
one can check that our setting recovers the standard multi-class logistic regression setting (see Foster et al. (2018) for an equivalent convention). We check below Assumptions (A1)–(A4).
Parameters λ, µ, β, ε, δ > 0,m ≥ 1 initialize A0 = λI, b0 = 0, θ1 = 0 for t = 1, ..., n do
receive xt ∈ Rd sample independently ω1, ..., ωm ∼ N (θ>t Φ(xt),Φ(xt)>A−1t−1Φ(xt)) predict ỹt = σ+(smoothµ( 1m ∑m i=1 σ(ωi))) where σ(z)i = ezi∑ j e zj , σ+(z)i = log(zi) and
smoothµ(z) := (1− µ)z + µ1/K. receive yt ∈ {1, ...,K} update θt+1 = argminθ∈RKd b > t−1θ + θ
>At−1θ + `xt,yt(θ) bt = bt−1 +∇`xt,yt(θt+1)− β∇2`xt,yt(θt+1)θt At = At−1 + β 2∇
2`xt,yt(θt+1) end
Algorithm 1: Efficient-GAF for K-class logistic regression
(A1) It holds with α = 1, since by Proposition 1 of Foster et al. (2018), the logistic loss is 1-mixable. Indeed, given a distribution π on RK , the choice ŷπ = σ+(Eŷ∼πσ(ŷ)), where σ+(z)k := log(zk), satisfies
Eŷ∼π exp(−`(ŷ, y)) = Eŷ∼πσ(ŷ)y = σ(ŷπ)y = exp(−`(ŷπ, y)) , for any y ∈ [K].
(A2) It is true for ζ = 4R2 by Lemma 6 in Appendix C. (A3) It holds with β = (log(K)/2 +BR+ 1)−1 by Lemma 4 in Appendix C applied with the
choices a = θ>1 Φ(x) and b = θ > 2 Φ(x) ∈ [−BR,BR]K .
(A4) It is valid with γ = R2. Indeed, ∇2`x,y(θ) = Φ(x)(diag(p) − pp>)Φ(x)> where p = σ(θ>Φ(x)) which can be bounded by ∇2`x,y(θ) ≤ Φ(x)Φ(x)> ≤ R2I . In particular, the proof of the mixability assumption (A1) provides us a concrete expression for Equation (4) of GAF by setting
ŷt = σ +(Eθ∼P̃t−1(θ)(σ(θ >xt))) . (6)
Substituting the above parameters into Theorem 1 yields the following corollary. Corollary 2. Let R,B > 0 and d′, n,K ≥ 1. Let (x1, y1), ..., (xn, yn) ∈ X × Y be an arbitrary sequence of observations and `(ŷ, y) = − log (σ(ŷ)y). GAF, run with λ = 32d′KR2, β = (log(K)/2 +BR+ 1)−1, and α = 1, satisfies
Rn(θ) . d′KB2R2 + d′K [log(K) +BR] log ( 1 + nd′K(log(K)+BR) ) , (7)
for all θ ∈ B(Rd, B), where . denotes an approximate inequality which is up to universal multiplicative and additive constants.
Computationally efficient approximation Yet, a key difficulty remains to compute exactly ŷt in Equation (6): the calculation of the expectation. In general, there is no closed-form expression and an approximation algorithm, such as Monte Carlo sampling, must be used. We provide now a fully implementable approximated version of GAF, which satisfies the same regret guarantees up to negligible additive constants. It is described in Algorithm 1.
First, to ensure that our approximation of the expectation produces forecasts close to those of the exact algorithm, we must smooth the function σ+ in Equation (6). Following the idea of Foster et al. (2018), we thus define for some µ ∈ [0, 12 ], the smoothing operator smoothµ : ∆K → ∆K by
smoothµ(p) = (1− µ)p+ µ1/K, ∀p ∈ ∆k . Now, we can show that σ+(smoothµ(·)) is µ-Lispchitz. Denoting the smoothed approximation of ŷt by ȳt = σ+(smoothµ(Eθ∼P̃t−1(σ(θxt))), Lemma 16 of Foster et al. (2018) shows that using ȳt instead of ŷt worsens the regret bound only by an additional constant µn; hence our choice µ = n−1. The last step of the approximated algorithm consists in using Monte Carlo sampling to approximate the expectation by a finite sum
ỹt = σ + ( smoothµ ( 1
m m∑ i=1 σ(ωi) )) where ωi i.i.d.∼ N (θ>t Φ(xt),Φ(xt)>A−1t−1Φ(xt)) .
Using Chernoff’s inequality, we show below that ỹt concentrates well to ȳt which entails the following guarantee.
Proposition 3. Let δ > 0. Efficient-GAF, run with parameters λ = 32d′KR2, β = (log(K)/2 + BR+ 1)−1, m ≥ 1, and µ = (log(n/δ)/m)1/3 satisfies with probability 1− δ the regret bound
Rn(θ) . d ′KB2R2 + d′K [log(K) +BR] log ( 1 +
n
d′K(log(K) +BR)
) +n ( K m log (n δ )) 13 with a computational cost O(nK3d2 +K2nm). In particular, the choice m = n3 log(n/δ), yields the regret bound (7) with probability 1 − δ and a total computational time of order O ( nK3d2 +
K2n4 log(n/δ) ) .
Proof. We first prove the computational complexity upper-bound. At each iteration t = 1, . . . , n, the algorithm performs the following computations:
(i) UpdateA−1t : since the rank ofAt−At−1 isK, the inverseA−1t can be updated inO(K3d2) operations.
(ii) Update θt+1: by Theorem 4 of Jézéquel et al. (2020), this can be done in O(log(n)+K2d2). (iii) Compute ỹt: to do so, it must sample m times from a Gaussian of dimension K. The cost is
thus O(mK2). Therefore, the overall time complexity of the algorithm is O ( nK3d2 +K2nm ) .
We now prove the corresponding regret bound. First, we define
p̃t := smoothµ ( 1 m ∑m i=1 σ(ωi) ) and p̄t := E[p̃t] = smoothµ ( Eθ∼P̃t−1(σ(θ >Φ(xt)) ) .
With this notation, ỹt = σ+(p̃t) and ȳt = σ+(p̄t). Using that σ(σ+(p)) = p for any p ∈ ∆K yields for all yt ∈ [K]
`(ỹt, yt) = − log ( σ(ỹt)yt ) = − log ( σ(σ+(p̃t))yt ) = − log(p̃t,yt) ,
where p̃t,i denotes the i-th component of p̃t. Similarly, `(ȳt, yt) = − log(p̄t,yt). Therefore, fixing some ε > 0, we have
P [ `(ỹt, yt)− `(ȳt, yt) > ε ] = P [ − log ( p̃t,yt p̄t,yt ) > ε ] = P [p̃t,yt − p̄t,yt < (e−ε − 1)p̄t,yt ] .
Using that e−ε − 1 ≤ 1− 2ε for ε ∈ (0, 12 ) together with the multiplicative Chernoff’s bound (see e.g., Theorem 1.10.5 of Doerr (2020)), entails for all 0 < ε < 12 ,
P [`(ỹt, yt)− `(ȳt, yt) > ε] ≤ P [p̃t,yt < (1− 2ε)p̄t,yt ] ≤ e−2ε 2mp̄t,yt ≤ e−ε 2mµ/K , (8)
where the last inequality is because smoothµ(p)i ≥ µ/K for all i ∈ [K]. Taking ε =√ K µm log(n/δ), we have
P [ `(ỹt, yt)− `(ȳt, yt) > √ K µm log ( n δ )] ≤ δn . (9)
Furthermore, denoting p̂t = Eθ∼P̃t−1(σ(θ >Φ(xt)) so that p̄t = smoothµ(p̂t) and following Lemma 16 of Foster et al. (2018), we get
`(ȳt, yt)− `(ŷt, yt) = log ( p̂t,yt p̄t,yt ) = log ( p̂t,yt smoothµ(p̂t)yt ) = log ( p̂t,yt (1− µ)p̂t,yt + µ K ) = log ( 1 1− µ ( 1− 1Kp̂t,yt )) ≤ 2µ(1− 1 Kp̂t,yt ) ≤ 2µ . (10)
Combining (9) and (10), and using a union bound yields with probability at least 1− δ∑n t=1 `(ỹt, yt)− `(ŷt, yt) ≤ n √ K µm log ( n δ ) + µn .
Optimizing µ = (log(n/δ)K/m)1/3 and using Corollary 3 concludes the proof.
Using Gaussian distribution improves significantly the computation time with respect to log-concave distribution as considered by Foster et al. (2018). Indeed, it allows to get exact and efficient samples from P̃t−1 which leads to a computation time of order O(n4). On the contrary, since it is not possible to draw an exact sample from a log-concave distribution, Foster et al. (2018) must resort to expensive random walks to approximately sample from it. Using the method from Bubeck et al. (2018), as suggested by Foster et al. (2018), leads to a O(B6n24(Bn+ dK)12) computation time per iteration (see their Example 3 with the choices ε = n−2 and µ = n−1).
Note that our analysis shows a clean trade-off between the computational time and the regret bound. One could set a smaller value for m at the price of a larger regret. In particular, it is worth pointing out that our analysis consider the worst case scenario. In the experiments that we considered, it seems that much less samples are sufficient to reach a good accuracy. This is partly because the lower-bound p̄t,yt ≥ µ ≈ 1/n in Inequality (8) may be very coarse for many t. If the p̄t,yt were of order Ω(1), one could choose m = O(n2) instead of m = O(n3). Another possible explanation comes from the fact that, if At is large enough, the variance of the samples is small and the convergence is much faster. In our experiments, the choice of m = 100 already provides a good approximation. We leave to future work the possibility to improve the complexity in favorable scenarios.
Comparison with Jézéquel et al. (2020) In the binary case, a more efficient algorithm called AIOLI has been introduced in Jézéquel et al. (2020). The latter achieves a similar regret bound as GAF, while having a computational complexity of only O(n(d2 + log(n))). This was possible because AIOLI does not rely on Monte Carlo sampling at all and uses only convex optimization instead. Therefore, one may wonder if it is possible to extend directly AIOLI to the multiclass setting avoiding Monte Carlo methods used in this work. In fact, we first tried to analyse the regret of natural extensions of AIOLI but we found the following inherent difficulties. On the intuitive side, AIOLI was based on the observation that when θt is far from 0 (let say θt 0) either yt = 1 and the curvature was advantageous or yt = −1 and θ̂t tends to the oracle θt+1. This intuition is a bit lost in the multiclass setting as several oracles are possible if yt 6= 1. On a more technical side, the analysis in Jézéquel et al. (2020) crucially relies on the relation g−ytt = −(1 +BR)ηtgt (Equation 20) which seems to have no equivalent in the multiclass setting. Thus, it remains an open question if an extension of AIOLI to the multiclass setting is possible.

heading: 5 Experiments
text: Although GAF is primarily theoretically motivated in a worst-case analysis, here we study its performance on real data sets. We consider three datasets (vehicle, shuttle, and segmentation taken from LIBSVM Data 1) and compare the performance of GAF with two well-used algorithms: Online Gradient Descent (OGD) (Zinkevich, 2003) and Online Newton Step (ONS) (Hazan et al., 2007). The algorithm of Foster et al. (2018) is not considered because of prohibitive computational complexity. Concerning the hyper-parameters, the values suggested by the theory are generally too conservative. We thus choose the best ones in a grid for each algorithm (λ, β ∈ [0.01, 0.03, 0.1, 0.3, 1, 3, 10]). The averaged losses over time are reported in Figure 1. We can remark that the performance of GAF is similar to the one of ONS when the number of samples is high. However, the learning of GAF seems more stable than ONS which leads sometimes (vehicle and segment) to better performance when there are few samples. This is not surprising as aggregating forecasters are hedging against worst-case scenario. Those results show that GAF is not only a theoretical oriented algorithm but could also be used successfully in practice. However, we expect it to perform best in a hard adversarial regime (close to the one described in Hazan et al. (2014)).
Conclusion and future work
We have shown how to leverage both mixability and quadratic approximations to design an algorithm GAF which achieves high statistical performance while being more efficient than existing algorithms. In particular, it achieves a new trade-off for multiclass logistic regression.
Some interesting questions are still remaining and left for future work. The linear dependence in B although better than what can be achieved by any proper algorithm is still sub-optimal compared
1https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/multiclass.html
to Foster et al. (2018) (logarithmic dependence). It may be possible to improve it by using other surrogates than quadratics. The essential point would be to prove an equivalent of Lemma 4 for those surrogates. The computational complexity of computing Equation (6) should also remain low. The computational complexity may also be improved. Like previous point, other surrogates could be used to make Equation (6) easier to compute. Finally, we believe that GAF should easily be extended to kernels. One would need to adapt the analysis to depend on the effective dimension (Jézéquel et al., 2019) of the RKHS instead of d. It would of interest to see how classical approximation algorithms, like Nystrom and random Fourier features, deteriorate statistical performance. Finally, we have only provided two examples of loss functions (multi-class logistic loss and squared loss). An intriguing question would be to see if the algorithm can be used to improve existing results for other losses such as Huber loss.

Review summaries:
The authors propose a new computationally efficient algorithm for multi-class logistic regression. The propose method is based on a tight quadratic approximation of the loss function (extended from Lemma 5 of Jezequel et al.) and utilizing the mixability of the logistic loss (and other properties) to compute the prediction like AA of vovk. They show that their regret bound does not have the infamous exponential constant while being computationally efficient (n^4). 

In this paper, the authors consider some online prediction/convex optimization problems, where it is assumed that the losses satisfy a mixability condition, as well as additional assumptions inspired by logistic regression. While these problems can be addressed by Vovk's aggregating algorithm, the resulting procedure would be computationally expensive and involved. Efficient and more direct procedures have been proposed for least-squares regression by Vovk and Azoury-Warmuth (VAW), and more recently for (binary) logistic regression by Mourtada and Gaiffas (MG) and Jezequel, Gaillard and Rudi (JGR). The present paper sets out to extend these constructions; the considered losses are assumed to be mixable, and satisfy conditions such as smoothness and a (weakened) self-concordance satisfied by logistic regression (lines 215-223). A procedure for this setting, called Gaussian aggregation forecaster (GAF), is proposed.  The main outcomes are the following: - In Sec 4.1, for linear regression the authors recover the VAW procedure and its guarantee. [However, there seems to be a possible issue with this derivation.] - In Sec 4.2, the authors consider multi-class logistic regression, and obtain a procedure and guarantee for this problem. These extend those of MG and JGR for binary logistic regression to the multi-class setting. They also go beyond what can be achieved by procedures such as online gradient descent, online Newton step or follow-the-regularized-leader. This is the main result of this work.

This paper considers the online convex optimization setting with mixable losses. The exponential weights algorithm is known to have constant regret with mixable losses, but it is also known to have prohibitively high running time for certain mixable losses. A mixable loss function of interest is the logistic loss used in logistic regression, for which faster algorithms such as Online Newton Step (ONS) suffer a double exponentially larger regret bound compared to exponential weights.  As a solution the authors propose a mix between ONS and exponential weights called GAF. To update, the algorithm computes a quadratic surrogate loss akin to the surrogate loss used in ONS, but uses the actual hessian rather than an approximation to the hessian to compute the surrogate loss. Another interesting trick that the authors use is that the point around which the quadratic surrogate loss is defined: it is the minimizer of the $l_2$ regularization plus the loss of the previous round and the surrogate losses of all rounds preceding the previous round.  The predictions of the algorithms are the same predictions that would be used for exponential weights. While for arbitrary priors this would pose a computational problem, the algorithm uses a gaussian prior. As is observed in prior work, with a gaussian prior and quadratic losses the posterior of exponential weights is also a gaussian. This allows the authors to devise more efficient predictions.   

With online learning for multiclass logistic regression, its regret has a trade-off against the computational cost. For example, ONS can be computed in $O(n)$ time ($n$: number of instances) while the regret is $O(e^B \log n)$: exponentially increased against the radius of regret computation $B$. [Foster 2018] can provide the regret $O(\log(Bn))$ but requires $O(n^{37})$ time. The proposed method can provide another trade-off between them, and achieved the regret $O(B^2 + B\log n)$ with the computational cost of $O(n^4)$. This is much faster than [Foster 2018] with acceptable computational cost.