Bayesian Model Averaging with Exponentiated Least Squares Loss

arXiv:1408.1234v3 [] 27 Feb 2018

Dong Dai
IMS Health, PA, USA

Lei Han
Tencent AI Lab, Shenzhen, China

Ting Yang
Rutgers University, NJ, USA

Tong Zhang
Tencent AI Lab, Shenzhen, China

Abstract
The model averaging problem is to average multiple models to achieve a prediction accuracy
not much worse than that of the best single model in terms of mean squared error. It is
known that if the models are misspecified, model averaging is superior to model selection.
Specifically, let n be the sample size, then the worst case regret of the former
√ decays at a rate
of O(1/n) while the worst case regret of the latter decays at a rate of O(1/ n). The recently
proposed Q-aggregation algorithm [1] solves the model averaging problem with the optimal
regret of O(1/n) both in expectation and in deviation; however it suffers from two limitations:
(1) for continuous dictionary, the proposed greedy algorithm for solving Q-aggregation is not
applicable; (2) the formulation of Q-aggregation appears ad hoc without clear intuition. This
paper examines a different approach to model averaging by considering a Bayes estimator for
deviation optimal model averaging by using exponentiated least squares loss. We establish a
primal-dual relationship of this estimator and that of Q-aggregation and propose new algorithms
that satisfactorily resolve the above mentioned limitations of Q-aggregation.

1

Introduction

This paper considers the model averaging problem, where the goal is to average multiple models
in order to achieve improved prediction accuracy. Let x1 , . . . , xn be n given design points from
a space X , let H = {f1 , . . . , fM } be a given dictionary of real valued functions on X and denote
f j = (fj (x1 ), . . . , fj (xn ))> ∈ Rn for each j. The goal is to estimate an unknown regression function
η : X → R at the design points based on observations
yi = η(xi ) + ξi ,
where ξ1 , . . . , ξn are i.i.d. variables from N (0, σ 2 ).
The performance of an estimator η̂ is measured by its mean squared error (MSE) defined by
n

MSE(η̂) =

1X
(η̂(xi ) − η(xi ))2 .
n
i=1

We want to find an estimator η̂ that mimics the function in the dictionary with the smallest
MSE. Formally, a good estimator η̂ should satisfy the following exact oracle inequality in a certain
probabilistic sense:
MSE(η̂) ≤ min MSE(fj ) + ∆(n, M, σ 2 ) ,
(1)
j=1,...,M

1

where the remainder term ∆ > 0 should be as small as possible.
The problem of model averaging has been well-studied, and it is known (see, e.g., [2, 3]) that the
smallest possible order for ∆(n, M, σ 2 ) is σ 2 log M/n for oracle inequalities in expectation, where
“the smallest possible” is understood in the following minimax sense. There exists a dictionary
H = {f1 , . . . , fM } such that the following lower bound holds. For any estimator η̂, there exists a
regression function η such that
E MSE(η̂) ≥

min MSE(fj ) + Cσ 2

j=1,...,M

log M
n

for some positive constant C. It also implies that the lower bound holds not only in expectation
but also with positive probability.
Although our goal is to achieve an MSE as close as that of the best model in H, it is known (see
Theorem 2.1 of [4]) that there exists a dictionary H such that any estimator η̂ taking values restricted to the elements of H (such an estimator is referred to as a model selection estimator)
cannot
p
achieve an oracle inequality of form (1) with a remainder term of order smaller than σ (log M )/n;
in other words, model selection is suboptimal for the purpose of competing with the best single
model from a given family.
Instead of model selection, we can employ model averaging to derive oracle inequalities of form
(1) that achieves the optimal regret in expectation (see the references in [4]). More recently,
several work has produced optimal oracle inequalities for model averaging that not only hold in
expectation but also in deviation [5, 6, 7, 8, 3, 1]. In particular, the current work is closely related to
the Q-aggregation estimator investigated in [1] which solves the optimal model averaging problem
both in expectation and in deviation with a remainder term ∆(n, M, σ 2 ) of order O(1/n); the
authors also proposed a greedy algorithm GMA-0 for Q-aggregation which improves the Greedy
Model Averaging (GMA) algorithm firstly proposed by [8]. Yet there are still two limitations of
Q-aggregation: (1) Q-aggregation can be generalized for continuous candidates dictionary H, but
the greedy model averaging method GMA-0 can not be applied in such setting; (2) Q-aggregation
can be regarded intuitively as regression with variance penalty, but it lacks a good theoretical
interpretation.
In this paper we introduce a novel method called Bayesian Model Averaging with Exponentiated
Least Squares Loss (BMAX). We note that the previously studied exponential weighted model aggregation estimator EWMA (e.g., [4]) is the Bayes estimator under the least squares loss (posterior
mean), which leads to optimal regret in expectation but is suboptimal in deviation. In contrast,
the new BMAX model averaging estimator is essentially a Bayes estimator under an appropriately
defined exponentiated least squares loss, and we will show that the Q-aggregation formulation (with
Kullback-Leibler entropy) in [1] is essentially a dual representation of the newly introduced BMAX
formulation, and it directly implies the optimality of the aggregate by BMAX. Computationally,
the new model aggregation method BMAX can be approximately solved by a greedy model averaging algorithm and a gradient descend algorithm which is applicable to continuous dictionary. In
summary, this paper establishes a natural Bayesian interpretation of Q-aggregation, and provides
additional computational procedures that are applicable for the continuous dictionary setting. This
relationship provides deeper understanding for model averaging procedures, and resolves the above
mentioned limitations of the Q-aggregation scheme.

2

2

Notations

This section introduces some notations used in this paper. In the following, we denote by Y =
(y1 , . . . , yn )> the observation vector, η = (η(x1 ), . . . , η(xn ))> the model output, and ξ = (ξ1 , . . . , ξn )>
the noise vector. The underlying statistical model can be expressed as
Y =η+ξ ,

(2)
P
with ξ ∼ N (0, σ 2 I n ). We also denote `2 norm as kY k2 = ( ni=1 yi2 )1/2 , and the inner product as
hξ, fi = ξ > f. Let ΛM be the flat simplex in RM defined by


M


X
ΛM = λ = (λ1 , . . . , λM )> ∈ RM : λj ≥ 0,
λj = 1 ,


j=1

and π = (π1 , . . . , πM )> ∈ ΛM be a given prior.
P
Each λ ∈ ΛM yields a model averaging estimator as fλ = M
j=1 λj fj ; that is, using the vector
PM
>
notation fλ = (fλ (x1 ), . . . , fλ (xn )) we have fλ = j=1 λj fj . The Kullback-Leibler divergence for
λ, π ∈ ΛM is defined as
K(λ, π) =

M
X

λj log(λj /πj ) ,

j=1

and in the definition we use the convention 0 · log(0) = 0. For matrices A, B ∈ Rn×n , A ≥ B
indicates that A − B is positive semi-definite.

3

Bayesian Model Averaging with Exponentiated Least Squares
Loss

The traditional Bayesian model averaging estimator is the exponential weighted model averaging
estimator EWMA [4] which optimizes the least squares loss. Although the estimator is optimal
in expectation, it is suboptimal in deviation [1]. In this section we introduce a different Bayesian
model averaging estimator called BMAX that optimizes an exponentiated least squares loss.
In order to introduce the BMAX estimator, we consider the following Bayesian framework,
where we should be noted that the assumptions below are only used to derive BMAX, and these
assumptions are not used in our theoretical analysis. Y is a normally distributed observation vector
with mean µ = (µ1 , . . . , µM )> and covariance matrix ω 2 I n :
Y |µ ∼ N (µ, ω 2 I n ) ,

(3)

and for j = 1, . . . , M , the prior for each model fj is
π(µ = fj ) = πj .
In this setting, the posterior distribution of µ given Y is


kfj −Y k22
πj
exp
−
2
p(Y |µ = fj )p(µ = fj )
2ω

 .
=P
p(µ = fj |Y ) = PM
kfj −Y k22
M
p(Y
|µ
=
f
)p(µ
=
f
)
j
j
exp
−
πj
j=1
j=1
2ω 2
3

(4)

In the Bayesian decision theoretical framework considered in this paper, the quantity of interest is
η = EY , and we consider a loss function L(ψ, µ) which we would like to minimize with respect to
the posterior distribution. The corresponding Bayes estimator ψ̂ minimizes the posterior expected
loss from µ as follows:
ψ̂ = argmin E [L(ψ, µ)|Y ] .
(5)
ψ∈Rn

It is worth pointing out that the above Bayesian framework is only used to obtain decision theoretically motivated model averaging estimators (Bayesian estimators have good theoretical properties
such as admissibility, etc). In particular we do not assume that the model itself is correctly specified. That is, in this paper we allow misspecified models, where the parameters µ and ω 2 are not
necessarily equal to the true mean η and the true variance σ 2 in (2), and η does not necessarily
belong to the dictionary {f1 , . . . , fM }.
The Bayesian estimator of (5) depends on the underlying loss function L(·, ·). For example,
under the standard least squares loss L(ψ, µ) = kψ − µk22 , the Bayes estimator is the posterior
mean, which leads to the Exponential Weighted Model Aggregation (EWMA) estimator [4] :


PM
kf −Y k2
exp − j2ω2 2 πj fj
j=1


ψ `2 (ω 2 ) = P
.
(6)
kfj −Y k22
M
exp
−
π
j
j=1
2ω 2
This estimator is optimal in expectation [9, 10], but suboptimal in deviation [1].
In this paper, we introduce the following exponentiated least squares loss motivated from the
exponential moment technique for proving large deviation tail bounds for sums of random variables:


1−ν
2
L(ψ, µ) = exp
kψ − µk2 ,
(7)
2ω 2
where the parameter ν ∈ (0, 1). It is easy to verify that the Bayes estimator defined by (5) with
the loss function defined in (7) can be written as
ψ X (ω 2 , ν) = argmin J(ψ) ,

(8)

ψ∈Rn

where
J(ψ) =

M
X
j=1



1−ν
1
2
2
πj exp − 2 kfj − Y k2 +
kψ − fj k2 .
2ω
2ω 2

(9)

The estimator ψ X (ω 2 , ν) will be referred to as the Bayesian model aggregation estimator with
exponentiated least squares loss (BMAX).
To minimize J(ψ), it is equivalent to minimize log J(ψ). Lemma 1 below shows the strong
convexity and smoothness (under some conditions) of log J(ψ).
Given ν ∈ (0, 1) and ω > 0, we define
A1 =

1−ν
;
ω2

(10)

moreover, if the `2 -norm of every fj is bounded by a constant L ∈ R:
kfj k2 ≤ L , ∀ j = 1, . . . , M ,
4

(11)

we define A2 and A3 as


1−ν
1−ν 2 2
A2 = 2 +
L ,
ω
ω2




1−ν
1−ν 2 4
2
A3 =
L +
L .
ω2
ω2
Lemma 1. For any ψ ∈ Rn , define the Hessian matrix of log J(ψ) as ∇2 log J(ψ) =
then we have
∇2 log J(ψ) ≥ A1 I n .

(12)
(13)
∂ 2 log J(ψ)
,
∂ψ∂ψ >

(14)

If {f1 , . . . , fM } satisfies condition (11), then
∇2 log J(ψ) ≤ A2 I n ,

(15)

where A1 and A2 are defined in (10) and (12).

4

Dual Representation and Q-aggregation

In this section, we will show that the Q-aggregation scheme of [1] with the standard Kullback-Leibler
entropy solves a dual representation of the BMAX formulation defined by (8).
Given Y and {f1 , . . . , fM }, Q-aggregation fλQ is defined as follows:
fλQ =

M
X

λQ
j fj ,

(16)

j=1
Q >
M such that
where λQ = (λQ
1 , . . . , λM ) ∈ Λ

λQ ∈ argmin Q(λ) ,

(17)

λ∈ΛM

Q(λ) = kfλ − Y k22 + ν

M
X

λj kfj − fλ k22 + 2ω 2 Kρ (λ, π) ,

(18)

j=1

for some ν ∈ (0, 1), where the ρ-entropy Kρ (λ, π) is defined as
Kρ (λ, π) =

M
X



ρ(λj )
πj


,

(19)

t log ρ(t) is convex .

(20)

λj log

j=1

where ρ is a real valued function on [0, 1] satisfying
ρ(t) ≥ t ,

When ρ(t) = t, Kρ (λ, π) becomes K(λ, π), i.e., the Kullback-Leibler entropy. When ρ(t) = 1,
P
M
Kρ (λ, π) = M
j=1 λj log(1/πj ), a linear entropy in Λ , and in particular the penalty Kρ (λ, π) in
(18) becomes a constant when π is a flat prior.

5

Now, to illustrate duality, we shall first introduce a function T : Rn → R as


M

 ν
X
ν
T (h) = −
kh − Y k22 − 2ω 2 log 
πj exp − 2 kfj − hk22  ,
1−ν
2ω

(21)

j=1

and denote the maximizer of T (h) as
ĥ = argmax T (h) .

(22)

h∈Rn

Define function S : ΛM × Rn → R as
M

S(λ, h) = −

X
ν
λj kfj − hk22 + 2ω 2 K(λ, π) .
kh − Y k22 + ν
1−ν

(23)

j=1

It is not difficult to verify that for ν ∈ (0, 1), S(λ, h) is convex in λ and concave in h. The following
duality lemma states the relationship between ĥ and fλQ .
Lemma 2. When ρ(t) = t, we have the following result
Q(λ) = maxn S(λ, h),

T (h) = min S(λ, h).
λ∈ΛM

h∈R

min Q(λ) = min maxn S(λ, h) = maxn min S(λ, h) = maxn T (h),

λ∈ΛM

λ∈ΛM h∈R

h∈R λ∈ΛM

h∈R

where the equality is achieved at (λQ , ĥ). Moreover, we have
n
o
(λQ , ĥ) = A ∩ B,
where A and B are two hyper-surfaces in ΛM × Rn defined as


1−ν
1
M
n
fλ ,
A = (λ, h) ∈ Λ × R : h = Y −
ν
ν
)
(

exp − 2ων 2 kfj − hk22 πj
M
n
.
B = (λ, h) ∈ Λ × R : λj = PM

ν
2
i=1 exp − 2ω 2 kfi − hk2 πi

(24)

Lemma 2 states that, (λQ , ĥ) is the only joint of hyper-surfaces A and B, and the only saddle
point of function S(λ, h) over space ΛM × Rn .
With T (h) defined as in (21), we can employ the transformation h = ν1 Y − 1−ν
ν ψ, and it is
easy to verify that
T (h) = −2ω 2 log (J(ψ)) ,

(25)

where J(ψ) is defined in (9). Since J(ψ) is strictly convex, T (h) is strictly concave and ĥ is unique.
It follows that maximizing T (h) is equivalent to minimizing J(ψ), and thus
1
1−ν
Y −
ψ X (ω 2 , ν) .
ν
ν
We can combine this representation with
ĥ =

ĥ =

1
1−ν
Y −
f Q
ν
ν λ

from Lemma 2 to obtain ψ X (ω 2 , ν) = fλQ . Therefore, we directly have the following relationship.
6

Theorem 1. When ρ(t) = t,
ψ X (ω 2 , ν) = fλQ ,
where ψ X (ω 2 , ν) is defined by (8) and (9), and fλQ is defined by (16),(17) and (18).
Theorem 1 states that, when ρ(t) = t, Kρ (λ, π) becomes the Kullback-Leibler entropy, and
Q-aggregation with the Kullback-Leibler entropy leads to an estimator fλQ that is essentially a
dual representation of the BMAX estimator ψ X (ω 2 , ν). It follows that, ψ X (ω 2 , ν) should share the
same optimality (both in expectation and deviation) as fλQ in solving the model averaging problem
(optimality of fλQ is shown in Theorem 3.1 of [1] with more general Kρ (λ, π), where ρ(t) only needs
to satisfy the condition (20)).
However, unlike the primal objective function J(ψ) which is defined on Rn , the dual objective
function Q(λ) is defined on RM . When M is large or infinity, the optimization of Q(λ) is non-trivial.
Although greedy algorithms are proposed in [1], they cannot handle the standard KL-divergence;
instead, they can only work with the linear entropy where ρ(t) = 1; it gives a larger penalty than the
standard KL-divergence (and thus worse resulting oracle inequality), and it cannot be generalized
to handle continuous dictionaries (because in such case the linear entropy with ρ(t) = 1 will always
be +∞). Therefore, the numerical greedy procedures of [1] converge to a solution with a worse
oracle bound than that of the solution for the primal formulation considered in this paper.
The following two corollaries (Corollary 1 and Corollary 2) are listed for illustration convenience.
They are directly derived from Theorem 1 and the optimality of fλQ in Theorem 3.1 of [1], so we
omit their proofs.
Corollary 1. Assume that ν ∈ (0, 1) and ω 2 ≥
inequality holds
kψ X (ω 2 , ν) − ηk22 ≤ ν

M
X

σ2
min(ν,1−ν) .

For any λ ∈ ΛM , the following oracle

λj kfj − ηk22 + (1 − ν) kfλ − ηk22 + 2ω 2 K(λ, πδ) ,

(26)

j=1

with probability at least 1 − δ. Moreover,
2

Ekψ X (ω , ν) −

ηk22

≤ν

M
X

λj kfj − ηk22 + (1 − ν) kfλ − ηk22 + 2ω 2 K(λ, π)

j=1

≤ kfλ −

ηk22

+ν

M
X

λj kfj − fλ k22 + 2ω 2 K(λ, π) .

(27)

j=1

P
2
Corollary 1 implies that when the term ν M
j=1 λj kfj −fλ k2 and the divergence term K(λ, π) are
small, ψ X (ω 2 , ν) can compete with an arbitrary fλ in the convex hull with any λ ∈ ΛM . Actually,
we can obtain an oracle inequality that competes with the best single model, which is the situation
that λ is at a vertex of the simplex ΛM :
Corollary 2. Under the assumptions of Corollary 1, ψ X (ω 2 , ν) satisfies



1
2
2
2
2
kψ X (ω , ν) − ηk2 ≤ min
kfj − ηk2 + 2ω log
,
j∈1,...,M
πj δ

7

(28)

with probability at least 1 − δ. Moreover,
2

Ekψ X (ω , ν) −

ηk22


≤

kfj −

min

j∈1,...,M

ηk22



2

+ 2ω log

1
πj


.

(29)

2

σ
It is also worth pointing out that the condition ω 2 ≥ min(ν,1−ν)
implies that ω 2 is at least greater
than 2σ 2 (when ν = 1/2), and intuitively this inflation of noise allows the Bayes estimator to handle
misspecification of the true mean η, which is not necessarily included in the dictionary H. Similar
observations were found by [11].
Finally we note that in the Bayesian framework stated in this section, when we change the
underlying loss function L(ψ, µ) from the standard least squares loss to the exponentiated least
squares loss (7), Bayes estimator changes from EWMA which is optimal only in expectation to
BMAX which is optimal both in expectation and in deviation. The difference is that the least
squares loss only controls the bias, while the exponentiated least squares loss controls both bias
and variance (as well as higher order moments) simultaneously. This can be seen by using Taylor
expansion


exp

1−ν
kψ − µk22
2ω 2



1−ν
=1+
kψ − µk22 + (1/2)
2ω 2



1−ν
kψ − µk22
2ω 2

2
+ ··· .

Since deviation bounds require us to control high order moments, the exponentiated least squares
loss is naturally suited for obtaining deviation bounds.
It is also natural to extend Corollary 1 from discrete candidates dictionary H = {f1 , . . . , fM }
to infinite dictionary (M = ∞) as well as continuous dictionary. For example, given a matrix
X ∈ Rn×d , we may consider a continuous dictionary parameterized by vector w as HΩ = {fw : fw =
Xw ∈ Rn } where Ω = {w : w ∈ Rd }. Then, we have the following result.
Corollary 3. Assume ν ∈ (0, 1) and for any distribution p(w) over the parameter space Ω, the
σ2
divergence term K(p, π) is finite for given prior π(w). Then, if ω 2 ≥ min(ν,1−ν)
, we have the oracle
inequality
kψ̂ −

Z

ηk22

≤ν

kfw −

ηk22 p(w)

2

Z
dw + (1 − ν)

Ω

+ 2ω 2 K(p, πδ)

fw p(w) dw − η
Ω

(30)

2

with probability at least 1 − δ. Moreover,
Ekψ̂ −

ηk22

Z
≤ν

kfw −

ηk22 p(w)

dw + (1 − ν)
2

Z
fw p(w) dw − η
Ω

fw −

+ν
2

2
2

Z

Z

p(w) dw + 2ω 2 K(p, π) ,

fw p(w) dw
Ω

Ω

+ 2ω 2 K(p, π)

fw p(w) dw − η
Ω

Ω

≤

2

Z

2

where
Z
ψ̂ = argmin
ψ∈Rn

exp
Ω

!
1
1−ν
2
2
− 2 kfw − Y k2 +
kfw − ψk2 π(w) dw.
2ω
2ω 2

8

(31)

When the distribution p on w ∈ Rd is concentrated around a single model, for example, p
is an uniform distribution on a small ball {w : kw − w0 k2 ≤ r} for some small r > 0, we have
R
R
1
r2
1
2
2
2
Ω fw p(w) dw = fw0 and the term n Ω kfw − fw0 k2 p(w)dw ≤ n kXkF is small if n kXkF is bounded
(k · kF is the matrix Frobenius norm). Although
a small r would lead to a larger divergence term
R
K(p, π) and there is a tradeoff between Ω kfw − fw0 k22 p(w)dw and K(p, π), by dividing the number
of observations n on both sides of Eq. (31), the term n1 K(p, π) could also be small as long as n
is sufficiently large. Therefore, Corollary 3 implies that ψ̂ can compete with a single model fw0 ,
similar to the case of discrete dictionary.

5

Algorithms to solve BMAX

In this section, we propose two algorithms, the Greedy Model Averaging (GMA-BMAX) algorithm
and the Gradient Descent (GD-BMAX) algorithm, to solve BMAX. The convergence rates of both
algorithms will be shown. Specifically, denote k as the number of iterations in the algorithms,
GMA-BMAX algorithm has a converge rate of O(1/k) and GD-BMAX algorithm converges with
a geometric rate of O(q k ) for some q ∈ (0, 1). Oracle inequalities will be derived for the estimators
in the iterations for both the algorithms.
Strong convexity of log J(ψ) directly implies that the minimizer ψ X (ω 2 , ν) is unique. Moreover,
it implies the following proposition which shows that an estimator that approximately minimizes
log J(ψ) satisfies an oracle inequality slightly worse than that of ψ X (ω 2 , ν) in Corollary 1. This
result suggests that we can employ appropriate numerical procedures to approximately solve (8),
and Corollary 1 implies an oracle inequality for such approximate solutions.
Proposition 1. Let ψ̂ be an -approximate minimizer of log J(ψ) for some  > 0 that log J(ψ̂) ≤
minψ log J(ψ) + . Then, we have
kψ̂ − ηk22 ≤ kψ X (ω 2 , ν) − ηk22 + 2

p

2/A1 kψ X (ω 2 , ν) − ηk2 +

2
.
A1

Next, we present the numerical algorithms to solve the BMAX problem.

5.1

Greedy Model Averaging Algorithm (GMA-BMAX)

The GMA-BMAX algorithm given in Algorithm 1 is a greedy algorithm that adds at most one
function from the dictionary H at each iteration. This feature is attractive as it outputs a ksparse solution that depends on at most k functions from the dictionary after k iterations. Similar
algorithms for model averaging have appeared in [8] and [1].
The following proposition follows from the standard analysis in [12, 13, 14]. It shows that the
estimator ψ (k) from Algorithm 1 converges to ψ X (ω 2 , ν).
Proposition 2. For ψ (k) as defined in Algorithm 1 (GMA-BMAX), if {f1 , . . . , fM } satisfies condition (11), then
8A3
log J(ψ (k) ) ≤ log J(ψ X (ω 2 , ν)) +
.
(32)
k+3
Proposition 2 states that, after running the GMA-BMAX algorithm for k steps to obtain ψ (k) ,
the corresponding objective value log J(ψ (k) ) converges to the optimal objective value log J(ψ X (ω 2 , ν))
9

Algorithm 1 Greedy Model Averaging Algorithm (GMA-BMAX)
Input: Noisy observation Y , dictionary H = {f1 , . . . , fM }, prior π ∈ ΛM , parameters ν, ω.
Output: Aggregate estimator ψ (k) .
Let ψ (0) = 0;
for k = 1, 2, . . . do
2
Set αk = k+1
;
(k)
J = argminj log J(ψ (k−1) + αk (fj − ψ (k−1) ));
ψ (k) = ψ (k−1) + αk (fJ (k) − ψ (k−1) );
end for
at a rate O(1/k). Combining this result with Proposition 1, we obtain the following oracle inequality,
which shows that the regret of the estimator ψ (k) after
√ running k steps of GMA-BMAX converges
2
to that of ψ X (ω , ν) in Corollary 1 at a rate O(1/ k).
2

σ
Proposition 3. Assume ν ∈ (0, 1) and if ω 2 ≥ min(ν,1−ν)
. Consider ψ (k) as in Algorithm 1
(GMA-BMAX). For any λ ∈ ΛM , the following oracle inequality holds

kψ (k) − ηk22 ≤ ν

M
X

λj kfj − ηk22 + (1 − ν) kfλ − ηk22 + 2ω 2 K(λ, πδ)

(33)

j=1

s
+2

16A3
16A3
kψ X (ω 2 , ν) − ηk2 +
A1 (k + 3)
A1 (k + 3)

with probability at least 1 − δ. Moreover,
Ekψ (k) − ηk22 ≤ ν

M
X

λj kfj − ηk22 + (1 − ν) kfλ − ηk22

(34)

j=1

s
+ 2ω 2 K(λ, π) + 2
From Proposition 3, if ω 2 ≥
kψ

(k)

−

σ2
min(ν,1−ν) ,

ηk22

≤ kfj −

16A3
16A3
Ekψ X (ω 2 , ν) − ηk2 +
.
A1 (k + 3)
A1 (k + 3)

for any j = 1, . . . , M we have
ηk22



2

+ 2ω log

1
πj δ



√
+ O(1/ k)



√
+ O(1/ k) .

with probability at least 1 − δ, and
Ekψ

(k)

−

ηk22

≤ kfj −

ηk22

2

+ 2ω log



1
πj

When k → ∞, ψ (k) achieves the optimal deviation bound. However, it does not imply optimal
deviation bound for ψ (k) with small k, while the greedy algorithms described in [8] (GMA) and [1]
(GMA-0) achieve optimal deviation bounds for small k when k ≥ 2. The advantage of GMA-BMAX
is that the resulting estimator ψ (k) competes with any fλ with λ ∈ ΛM under the KL entropy,
and such a result can be applied even for infinity dictionaries containing functions indexed by
10

Algorithm 2 Gradient Descent Algorithm (GD-BMAX)
Input: Noisy observation Y , dictionary H = {f1 , . . . , fM }, prior π ∈ ΛM , parameters ν, ω 2 .
Output: Aggregate estimator ψ (k) .
Let ψ (0) = 0.
for k = 1, 2, . . . do
Choose step size tk ∈ (0, 2/A2 ) for k > 0;
Let
fλ(k−1) =

M
X

(k−1)

λj

fj ,

j=1

where λ(k−1) ∈ ΛM and
(k−1)
λj

∝ πj exp

!
1 − ν (k−1)
1
2
2
kψ
− fj k2 ;
− 2 kfj − Y k2 +
2ω
2ω 2

(35)

ψ (k) = (1 − tk 1−ν
)ψ (k−1) + tk 1−ν
f
;
ω2
ω 2 λ(k−1)
end for
continuous parameters, as long as the KL divergence K(λ, π) is well-defined (see relevant discussions
in Section 4). On the other hand, the greedy estimators of [1] for the Q-aggregation scheme can
only deal with an upper bound of KL divergence referred to as linear entropy (see Section 4) that
is not well-defined for continuous dictionaries. This means that GMA-BMAX is more generally
applicable than the corresponding greedy algorithm GMA-0 in [1].

5.2

Gradient Descent Algorithm (GD-BMAX)

An alternative way solving the BMAX is to use the gradient descend method. The GD-BMAX
algorithm is shown in Algorithm 2. The gradient is
∇ log J(ψ (k−1) ) =

∇J(ψ (k−1) )
J(ψ

(k−1)

)

=

1 − ν (k−1)
(ψ
− fλ(k−1) ) ,
ω2

where λ(k−1) ∈ ΛM is defined as (35). In the k-th step, the update operation is
ψ (k) = (1 − tk

1 − ν (k−1)
1−ν
)ψ
+ tk 2 fλ(k−1) = ψ (k−1) − tk ∇ log J(ψ (k−1) ) .
2
ω
ω

Therefore, Algorithm 2 is a gradient decent algorithm with step size tk . The following proposition
shows the convergence of the GD-BMAX algorithm.
Proposition 4. For ψ (k) as defined in Algorithm 2, if we choose a fixed step size tk = s ∈ (0, 2/A2 )
for k > 0 and {f1 , . . . , fM } satisfy condition (11), then


log J(ψ (k) ) − log J(ψ X (ω 2 , ν)) ≤ [1 − 2A1 (s − (A2 /2)s2 )]k log J(ψ (0) ) − log J(ψ X (ω 2 , ν)) .
(36)

11

Remark 1. We may choose tk = s = 1/A2 to minimize the righthand side of (36) such that


log J(ψ (k) ) − log J(ψ X (ω 2 , ν)) ≤ (1 − A1 /A2 )k log J(ψ (0) ) − log J(ψ X (ω 2 , ν)) .
Proposition 4 shows that the GD-BMAX algorithm converges at a geometric rate of O(q k ) with
q = 1 − 2A1 (s − (A2 /2)s2 ). Moreover, we have the following oracle inequality, which shows that the
regret of the estimator ψ (k) after running k steps of GD-BMAX converges to that of ψ X (ω 2 , ν) in
Corollary 1 at a rate of O(q k ).
Proposition 5. Assume ν ∈ (0, 1) and if ω 2 ≥
λ ∈ ΛM ,
kψ

(k)

−

ηk22

≤ν

M
X

σ2
min(ν,1−ν) ,

we have oracle inequality for any

λj kfj − ηk22 + (1 − ν) kfλ − ηk22 + 2ω 2 K(λ, π) + 2ω 2 log(1/δ)

(37)

j=1

+2

q
L2 [1 − 2A1 (s − (A2 /2)s2 )]k kψ X (ω 2 , ν) − ηk2 + L2 [1 − 2A1 (s − (A2 /2)s2 )]k ,

with probability at least 1 − δ. Moreover,
Ekψ (k) − ηk22 ≤ ν

M
X

λj kfj − ηk22 + (1 − ν) kfλ − ηk22 + 2ω 2 K(λ, π)

(38)

j=1

+2

q
L2 [1 − 2A1 (s − (A2 /2)s2 )]k Ekψ X (ω 2 , ν) − ηk2 + L2 [1 − 2A1 (s − (A2 /2)s2 )]k .

From Proposition 5, if ω 2 ≥
kψ

σ2
min(ν,1−ν) ,

(k)

−

ηk22

for any j = 1, . . . , M we have

≤ kfj −

ηk22

2

1
πj δ









+ 2ω log

+ O(q k ) ,

with probability at least 1 − δ and
Ekψ (k) − ηk22 ≤ kfj − ηk22 + 2ω 2 log

1
πj

+ O(q k ) ,

for some constant q = 1 − 2A1 (s − (A2 /2)s2 ) ∈ (0, 1).
Though the GD-BMAX algorithm does not give sparse output as the GMA-BMAX algorithm
does, it has a faster geometric convergence compared to GMA-BMAX. Similar to Proposition 3,
the results in Proposition 5 do not imply optimal deviation bounds of ψ (k) for small k (k < ∞).
The GD-BMAX algorithm can be naturally applied to continuous dictionary case. When the
dictionary is continuous or M is large, direct calculation of λ(k−1) ∈ ΛM in (35) is impractical
and we can use the Markov Chain Monte Carlo (MCMC) sampling method [15] to approximate
λ(k−1) for the k-th iteration in Algorithm 2. A Metropolis-Hastings (MH) algorithm is given in
(k−1)
Algorithm 3. The MH algorithm approximates fλ(k−1) with uT
, so the resulting estimator ψ (k)
at the k-th step in Algorithm 2 will have perturbations. Below we provide a result showing how
the perturbations from approximating fλ(k−1) would influence the convergence of log J(ψ (k) ) to
log J(ψ X (ω 2 , ν)).
12

Algorithm 3 Metropolis-Hastings (MH) sampler for estimating fλ(k−1) at the k-th step of Algorithm 2
Input: Noisy observation Y , dictionary H = {f1 , . . . , fM }, prior π ∈ ΛM , parameters ν, ω 2 ,
(k − 1)-th step estimator ψ (k−1) .
P
(k−1)
(k−1)
Output: uT
as estimator of fλ(k−1) = M
fj .
j=1 λj
Initialize j(0) = 0;
for t = 1, · · · , T0 + T do
Generate j̃ ∼ q(·|j (t−1) ) (e.g., q can be chosen as a Gaussian distribution with mean fj (t−1) );
Compute
!
q(j (t−1) |j̃)θ(j̃)
(t−1)
ρ(j
, j̃) = min
,1 ,
q(j̃|j (t−1) )θ(j (t−1) )
where
θ(j) = πj exp

!
1
1 − ν (k−1)
2
2
− 2 kfj − Y k2 +
kψ
− fj k2 ;
2ω
2ω 2

Generate a random variable u ∈ [0, 1] and let
(
j̃ ,
if u ≤ ρ(j (t−1) , j̃);
(t)
j =
j (t−1) , otherwise;
end for
Calculate
(k−1)

uT

=

1
T

TX
0 +T

fj (t) ;

t=T0 +1

(k−1)

Proposition 6. Given Y ∈ Rn , for all k > 0, we assume that uT
(k−1)

E[uT

kCOV

|ψ (k−1) ] = fλ(k−1) ,

(k−1)
[uT
|ψ (k−1) ]kop

≤ s2 ,

in Algorithm 3 satisfies
(39)
(40)

where k · kop is the matrix spectral norm. Then, we have




E log J(ψ (k) ) − log J(ψ X (ω 2 , ν)) ≤ [1 − 2A1 (s − (A2 /2)s2 )]k log J(ψ (0) ) − log J(ψ X (ω 2 , ν))
+ A1 ns2 /2 .
A variant of the GD-BMAX algorithm with continuous dictionary notation is also provided in
Algorithm 4, in which we assume the dictionary is parameterized by w as HΩ = {f(w) : f(w) ∈ Rn }
and Ω = {w : w ∈ Rd }.

6

Experiments

Although the contribution of this work is mainly theoretical, we include some simulations to illustrate the performance of the GMA-BMAX algorithm and GD-BMAX algorithm proposed for the
13

Algorithm 4 Gradient Descent Algorithm with Continuous Dictionary
Input: Noisy observation Y , continuous dictionary HΩ , prior π Ω , parameters ν, ω 2 .
Output: Aggregate estimator ψ (k) .
Let ψ (0) = 0.
for k = 1, 2, . . . do
Choose step size tk ∈ (0, 2/A2 ) for k > 0;
Initialize w(0) ;
for t = 1, · · · , T0 + T do
Generate w̃ ∼ q(·|w(t−1) );
Compute
γ(w

(t−1)

!
q(w(t−1) |w̃)θ(w̃)
,1 ,
q(w̃|w(t−1) )θ(w(t−1) )

, w̃) = min

where
θ(w) = π(w) exp

!
1
1 − ν (k−1)
2
2
− 2 kf(w) − Y k2 +
kψ
− f(w)k2 ;
2ω
2ω 2

Generate a random variable u ∈ [0, 1] and let
(
w̃ ,
if u ≤ γ(w(t−1) , w̃);
(t)
w =
w(t−1) , otherwise;
end for
Calculate
(k−1)

uT

ψ (k) = (1 − tk

=

1
T

TX
0 +T

f(w(t) );

t=T0 +1

1 − ν (k−1)
1 − ν (k−1)
)ψ
+ tk 2 uT
;
ω2
ω

end for
BMAX method. We focus on the average performance of different algorithms and configurations.
The simulations will focus on discrete dictionary in order to compare the BMAX method with
existing algorithms while the BMAX method can deal with continuous dictionary. Set n = 50 and
M = 500. We identify a function f with a vector (f (x1 ), . . . , f (xn ))> ∈ Rn . Let I n denote the
identity matrix of Rn and let Θ ∼ N (0, I n ) be a random vector, and define {f1 , . . . , fM } as
(
fj = Θ + s · ζ j for 1 ≤ j ≤ M1 ,
(41)
fj = ζ j for M1 < j ≤ M ,
where ζ j ∼ N (0, I n ) (j = 1, . . . , M ) are independent random vectors.
Let ∆ ∼ N (0, I n ) be a random vector. The regression function is defined by η = f1 + 0.5∆.
Note that typically f1 will be the closest function to η but not necessarily. The noise vector
ξ ∼ N (0, σ 2 I n ) is independent of {f1 , . . . , fM } and σ = 2.
14

We define the oracle model (OM) fk∗ , where k ∗ = argminj MSE(fj ). The model fk∗ is clearly
not a valid estimator because it depends on the unobserved η, however it can be used as a performance benchmark. The performance difference between an estimator η̂ and the oracle model fk∗
is measured by the regret defined as:
R(η̂) = MSE(η̂) − MSE(fk∗ ) .

(42)

Since the target is η = f1 + 0.5∆, and f1 and ∆ are random Gaussian vectors, the oracle model
is likely f1 (but it may not be f1 due to the misspecification vector ∆). The noise σ = 2 is relatively
large, which implies a situation where the best convex aggregation does not outperform the oracle
model. This is the scenario we considered here. For simplicity, all algorithms use a flat prior
πj = 1/M for all j. The experiment is performed with 100 replications.
One method compared is the STAR algorithm of [5], which is optimal both in expectation and
in deviation under the uniform prior. Mathematically, suppose fk1 is the empirical risk minimizer
among functions in H, where
[ j) ,
k1 = argmin MSE(f
(43)
j

the STAR estimator f ∗ is defined as
f ∗ = (1 − α∗ )fk1 + α∗ fk2 ,

(44)


[ (1 − α)fk1 + αfj , α ∈ (0, 1).
(α∗ , k2 ) = argmin MSE

(45)

where
α,j

Another natural solution to solve the model averaging problem is to take the vector of weights
λproj defined by
[ λ) ,
λproj ∈ argmin MSE(f
(46)
λ∈ΛM

which minimizes the empirical risk. We call λproj the vector of projection weights since the aggregate estimator fλproj is the projection of Y onto the convex hull of the fj ’s.
We compare the exponential weighted model averaging method denoted as EWMA. Q-aggregation
with linear entropy (when ρ(t) = 1) is also compared and will be solved by GMA-0 from [1] (see
Algorithm 5 below).
Algorithm 5 GMA-0 Algorithm
Input: Noisy observation Y , dictionary H = {f1 , . . . , fM }, prior π ∈ ΛM , parameters ν, β.
Output: Aggregate estimator fλ(k) .
Let λ(0) = 0, fλ(0) = 0;
for k = 1, 2, . . . do
2
;
Set αk = k+1
(k)
J = argminj Q(λ(k−1) + αk (e(j) − λ(k−1) ));
λ(k) = λ(k−1) + αk (e(J
end for

(k) )

− λ(k−1) );

We adopt flat priors π = 1/M (j = 1, . . . , M ) for simplicity. From the definition of Q(λ) (18),
it is easy to see that, the minimizer of Q(λ) (when ρ(t) = 1 with flat prior) becomes λproj in (46)

15

by setting ν = 0, so λproj is approximated by GMA-0 with ν = 0 and 200 iterations, and the
projection algorithm is denoted by PROJ.
We evaluate two versions of the GD-BMAX algorithm that one is exactly the algorithm depicted
in Algorithm 2 and the other is to approximate fλ(k−1) in Algorithm 2 with the MH sampler in
Algorithm 3. We denote the latter variant as the GD-MH-BMAX algorithm. For the MH sampler,
we use Gaussian distribution for q(·|·) and set T0 = T = 500. The GMA-BMAX, GD-BMAX
GD-MH-BMAX and GMA-0 algorithms are run for K iterations up to K = 150, with ν = 1/2.
Parameter ω for GMA-BMAX, GD-BMAX, GD-MH-BMAX and EWMA is tuned by 10-fold crossvalidation. Regrets of all algorithms defined in (42) are reported for comparisons.
In the following, we consider two scenarios. The first situation is when the bases are not very
correlated; in such case GMA-0 can perform better than GMA-BMAX, GD-BMAX and GD-MHBMAX because the former (which employs linear entropy) produces sparser estimators. The second
situation is when the bases are highly correlated; in such case GMA-BMAX, GD-BMAX and GDMH-BMAX are superior than GMA-0 because the former algorithms (which employ strongly convex
KL-entropy) give the clustered basis functions similar weights while the GMA-0 tends to select one
from the clustered basis functions, which may lead to model selection error. The correlated bases
situation occurs in the continuous dictionary setting.

6.1

Experiment 1: when s = 1 and M1 = 50, bases are not very correlated
Table 1: Performance Comparison (s = 1 and M1 = 50)
STAR
0.3895±0.3816

GMA-BMAX
GD-BMAX
GD-MH-BMAX
GMA-0

k=1
0.5726±0.7536
1.6187±0.3313
1.6192±0.3311
0.3742±0.7804

PROJ
0.3953±0.2680

k=5
0.4196±0.4339
0.7231±0.2305
0.7253±0.2328
0.2916±0.3691

k = 15
0.3480±0.3606
0.3003±0.3098
0.3061±0.3127
0.2567±0.3291

EWMA
0.2823±0.5718
k = 60
0.2840±0.3213
0.2603±0.3161
0.2645±0.3107
0.2546±0.3289

k = 100
0.2745±0.3180
0.2602±0.3161
0.2705±0.3237
0.2548±0.3292

k = 150
0.2690±0.3175
0.2602±0.3161
0.2583±0.3058
0.2550±0.3287

Table 2: Cumulative Frequency of Regret (s = 1, M1 = 50 and k = 150)
Upper Boundary
GMA-BMAX
GD-BMAX
GD-MH-BMAX
GMA-0
EWMA

0.00
17
13
14
21
39

0.33
69
70
68
72
78

0.66
85
85
85
84
82

0.98
95
96
97
96
84

1.31
100
100
100
100
88

1.64
100
100
100
100
95

1.97
100
100
100
100
99

2.29
100
100
100
100
100

Table 1 compares the commonly used estimators, i.e., STAR, PROJ and EWMA, with GMABMAX, GD-BMAX, GD-MH-BMAX, and GMA-0. The regrets are reported using the format of
“mean ± standard deviation”. Table 2 reports the cumulative frequency of GMA-BMAX, GDBMAX, GD-MH-BMAX, GMA-0 and EWMA with 100 replicates and fixed iteration k. For each
entry, we summarize the number of replicates with regrets which are smaller than or equal to the
upper boundary value.
16

The results in Table 1 indicate that GMA-0 achieves the best performance as iteration k increases. GMA-0 outperforms STAR, EWMA and PROJ after as small as k = 15 iterations, which
still gives a relatively sparse averaged model. This is consistent with Theorems 4.1 and 4.2 in [1]
which shows that GMA-0 has optimal bounds for small k (k ≥ 2).
GMA-BMAX, GD-BMAX and GD-MH-BMAX prefer dense model that assigns similar weights
to similar candidates. It is easy to verify that ψ X (ω 2 , ν) = fλ with λ ∈ ΛM defined as


1
1−ν
2
2
2
λj ∝ πj exp − 2 kfj − Y k2 +
kψ X (ω , ν) − fj k2 ,
2ω
2ω 2
and thus two similar candidates fi and fj will have similar weights.
In contrast, GMA-0 prefers sparse model by selecting candidates less related to the estimator
from previous iteration. Specifically, with flat prior π, the choice of J (k) in GMA-0 algorithm can
be further simplified to

J (k) = argmin kfj − Y k22 − (1 − ν)(1 − αk )kfλ(k−1) − fj k22 ,
(47)
j

and thus at each iteration k in the GMA-0 algorithm, the estimator fj is preferred if it is close to
Y while it is less correlated to the current aggregate estimator fλ(k−1) (because the minimization
requires kfλ(k−1) − fj k22 to be large while kfj − Y k22 being small).
In Experiment 1, the first 50 candidates {f1 , . . . , fM } are closer to the truth η = f1 + 0.5∆
than other candidate fj for j > 50, yet they are not very correlated when s = 1. Sparsity is
preferred when correlations are not strong among the predicting features (in our experiment, the
first 50 candidates), and GMA-0 is to output sparser estimator than GMA-BMAX, GD-BMAX and
GD-MH-BMAX. Therefore, we would expect GMA-0 achieving smaller regret than GMA-BMAX,
GD-BMAX and GD-MH-BMAX under this situation. Although GMA-BMAX, GD-BMAX and
GD-MH-BMAX are worse than GMA-0 when the bases are not very correlated, they beat EWMA,
STAR and PROJ when the iteration k is large enough.
Figure 1 compares the regrets of GMA-BMAX, GD-BMAX, GD-MH-BMAX, GMA-0 and
EWMA. (a-e) illustrate the histograms of the regrets with 100 replicates. The corresponding cumulative frequencies are presented in Table 2. Since Proposition 3 indicates that the optimal deviation
bound is obtained by k → ∞, we pick k = 150 for GMA-BMAX, GD-BMAX, GD-MH-BMAX and
GMA-0 in order to make fair comparison. As the histograms show, although EWMA has the most
replicates in which the regrets are close to zero, the distribution of the EWMA estimator is the most
dispersive with many extreme values compared to other methods. The performance is consistent
with [9, 10] and [1] which state that the EWMA estimator is optimal in expectation but sub-optimal
in deviation. Therefore, we would expect GMA-BMAX, GD-BMAX, GD-MH-BMAX and GMA-0
enjoy more concentrated distribution than EWMA because they are also optimal in deviation. (f)
shows the convergence of GMA-BMAX, GD-BMAX, GD-MH-BMAX and GMA-0. We observe
that GD-BMAX and GD-MH-BMAX show faster convergence compared to GMA-BMAX, which is
consistent with Proposition 5. Moreover, the GD-MH-BMAX algorithm approximates GD-BMAX
well with small perturbations, which is consistent with Proposition 6.
Note that GMA-BMAX and GMA-0 produce different estimators after the first iteration (k = 1).
GMA-BMAX selects j ∈ {1, . . . , M } that minimizes log J(fj ), while GMA-0 selects j ∈ {1, . . . , M }
that minimizes Q(fj ) and the output of the first stage is actually the empirical risk minimizer fk1
[ j ). Moreover, since sparse model is preferred in this scenario, GMA-0
where k1 = argminj MSE(f
has smaller regret than GMA-BMAX, GD-BMAX and GD-MH-BMAX, and all of them converge
fast within a few iterations.
17

GMA-BMAX

GD-BMAX

80

70

60

60

60

50

50

50

30

Frequency

70

40

40

40

30

30

20

20

20

10

10

10

0

0

0.5

1

1.5

0

2

0

0

0.5

1

Regret

1.5

2

0

0.5

1

(b)

GMA-0

70

60

60

50

50

2

(c)

EWMA

80

70

1.5

Regret

Regret

(a)
80

GD-MH-BMAX

80

70

Frequency

Frequency

80

1.6
GMA-BMAX
GD-BMAX
GD-MH-BMAX
GMA-0

1.4

40
30

Regret

Frequency

Frequency

1.2

40

1
0.8

30

20

20

0.6

10

10

0.4

0

0

0

0.5

1

1.5

2

0

0.5

Regret

1

1.5

2

0

(d)

50

100

150

Iteration

Regret

(e)

(f)

Figure 1: (a-e) show the histograms of regrets for GMA-BMAX, GD-BMAX, GD-MH-BMAX,
GMA-0 (with k = 150) and EWMA; (f) reports the results of regrets R(ψ (k) ) versus iterations k
under the case s = 1 and M1 = 50.

6.2

Experiment 2: when s = σ/kζ j k and M1 = M , bases are all highly correlated
Table 3: Performance Comparison(s = σ/kζ j k and M1 = M )
STAR
0.0587±0.0338

GMA-BMAX
GD-BMAX
GD-MH-BMAX
GMA-0

k=1
0.0929±0.0381
0.7817±0.2207
0.7414±0.2126
0.0904±0.0380

k=5
0.0503±0.0300
0.2525±0.1041
0.1996±0.0876
0.0650±0.0379

PROJ
0.0495±0.0322
k = 15
0.0431±0.0285
0.0501±0.0317
0.0435±0.0291
0.0630±0.0371

EWMA
0.0452±0.0307
k = 60
0.0408±0.0282
0.0406±0.0284
0.0393±0.0273
0.0621±0.0364

k = 100
0.0407±0.0283
0.0406±0.0284
0.0382±0.0273
0.0622±0.0364

k = 150
0.0406±0.0283
0.0406±0.0284
0.0384±0.0284
0.0620±0.0364

In Experiment 2, we define regression function as η = Θ + 0.5∆ which is slightly different
from Experiment 1. The results in Table 3 indicate that GMA-BMAX, GD-BMAX and GD-MHBMAX perform better than GMA-0 as iteration k increases, and GMA-BMAX, GD-BMAX and
GD-MH-BMAX also beat STAR, PROJ and EWMA when k is large enough.
In this experiment, all of the candidates {f1 , . . . , fM } are close to the truth η = Θ + 0.5∆, and
18

Table 4: Cumulative Frequency of Regret (s = σ/kζ j k, M1 = M and k = 150)
Upper Boundary
GMA-BMAX
GD-BMAX
GD-MH-BMAX
GMA-0
EWMA

0.000
5
6
5
4
5

0.026
31
31
36
17
26

0.051
69
69
71
40
59

0.077
89
88
90
67
86

0.103
97
97
97
89
96

0.128
99
99
99
95
98

0.154
100
100
100
99
100

0.180
100
100
100
100
100

they are highly correlated when s is small. GMA-BMAX, GD-BMAX and GD-MH-BMAX tend
to assign similar weights to similar candidates, while GMA-0 tends to exclude other correlated
candidates once one has been selected. Therefore, GMA-BMAX, GD-BMAX and GD-MH-BMAX
will average over those candidates with similar weights resulting less variance (also less bias due
to the design), while GMA-0 will have high variance because it only selects one of the candidates.
Moreover, similar to GMA-BMAX, GD-BMAX and GD-MH-BMAX, EWMA also prefer dense
model by assigning similar weights to similar candidates. However, the EWMA estimator assigns
weights defined as


1
λj ∝ πj exp − 2 kfj − Y k22 ,
2ω
while the weights in GMA-BMAX, GD-BMAX and GD-MH-BMAX are defined as


1
1−ν
2
2
λj ∝ πj exp − 2 kfj − Y k22 +
kψ
(ω
,
ν)
−
f
k
,
j 2
X
2ω
2ω 2
where ψ X (ω 2 , ν) = fλ with λ ∈ ΛM . Note that for all j, kfj − Y k22 are roughly equal to each
other under this scenario. That is, the EWMA estimator becomes the average of all candidates.
However, the GMA-BMAX, GD-BMAX and GD-MH-BMAX estimators are still weighed averages
of all bases, and the weights are adjusted by the extra term kψ X (ω 2 , ν) − fj k22 . Therefore, we would
hope GMA-BMAX, GD-BMAX and GD-MH-BMAX have smaller variances than EWMA.
Figure 2 compares the regrets of GMA-BMAX, GD-BMAX, GD-MH-BMAX, GMA-0 and
EWMA. (a-e) summarize the histograms of the regrets, and the corresponding cumulative frequencies are represented in Table 4. GMA-BMAX, GD-BMAX and GD-MH-BMAX have the most
concentrated distributions, because they are optimal both in expectation and in deviation. (f)
illustrates the convergence of GMA-BMAX, GD-BMAX, GD-MH-BMAX and GMA-0. All these
methods converge within a few iterations. As expected, GMA-BMAX, GD-BMAX and GD-MHBMAX achieve lower regret than GMA-0 when the basis functions are correlated.

7

Conclusion

This paper introduces a new formulation for deviation optimal model averaging which we refer to
as BMAX. It is motivated by Bayesian theoretical considerations with an appropriately defined exponentiated least squares loss. Moreover we established a primal-dual relationship of this estimator
and the Q-aggregation scheme (with KL entropy) by [1]. This relationship not only establishes a
natural Bayesian interpretation for Q-aggregation but also leads to new numerical algorithms for

19

GMA-BMAX

GD-BMAX

20

18

18

16

16

18

14

14

12

12

8

16
14

Frequency

10

10
8

12
10
8

6

6

4

4

4

2

2

2

0

0

0

0.05

0.1

0.15

6

0

0

0.05

Regret

0.1

0.15

0

0.05

(b)

GMA-0

18

16

16

14

14

12

12

0.15

(c)

EWMA

20

18

0.1

Regret

Regret

(a)
20

GD-MH-BMAX

22
20

Frequency

Frequency

20

GMA-BMAX
GD-BMAX
GD-MH-BMAX
GMA-0

0.7

10
8

0.5

Regret

Frequency

Frequency

0.6

10

0.4

8
0.3

6

6

4

4

0.2

2

2

0.1

0

0

0

0.05

0.1

0.15

0

0.05

Regret

0.1

0

0.15

50

(d)

100

150

Iteration

Regret

(e)

(f)

Figure 2: (a-e) show the histograms of regrets for GMA-BMAX, GD-BMAX, GD-MH-BMAX,
GMA-0 (with k = 150) and EWMA; (f) reports the result of regrets R(ψ (k) ) versus iterations k
under the case s = σ/kζ j k and M1 = M .
model aggregation that are suitable for the continuous dictionary setting where some basis functions are highly correlated. The new formulation and its relationship to Q-aggregation provides
deeper understanding of deviation optimal model averaging procedures.

A

Proofs

A.1

Proof of Lemma 1

Define λ ∈ ΛM as



1
1−ν
2
2
λj ∝ πj exp − 2 kfj − Y k2 +
kψ − fj k2 .
2ω
2ω 2
It follows that
∇J(ψ)
1−ν
=
(ψ − fλ ),
J(ψ)
ω2
and
M

∇2 J(ψ) X
=
λj
J(ψ)
j=1



1−ν
ω2

2

>

(ψ − fj )(ψ − fj ) +

20



1−ν
ω2

!


In

.

Then we have
(∇2 J(ψ))J(ψ) − (∇J(ψ))(∇J(ψ))>
J 2 (ψ)


 ! 


M
X
1−ν 2
1−ν 2
1
−
ν
>
=
In −
λj
(ψ − fj )(ψ − fj ) +
(ψ − fλ )(ψ − fλ )>
ω2
ω2
ω2

∇2 log J(ψ) =

j=1


=

1−ν
ω2


In +

M
X


λj

j=1

1−ν
ω2

2

(fλ − fj )(fλ − fj )> .


Therefore, ∇2 log J(ψ) ≥ 1−ν
I n.
ω2
With the assumption that kfj k2 ≤ L for all j, we have
M
X

λj (fλ − fj )(fλ − fj )> =

j=1

λj fj fj> − fλ fλ> ≤

j=1

It follows that ∇2 log J(ψ) ≤

A.2

M
X



1−ν
ω2



+

M
X

λj fj fj> ≤

j=1


1−ν 2 2
L
ω2



M
X

λj L2 I n = L2 I n .

j=1

I n.

Proof of Lemma 2

Lemma 3. For any λ ∈ ΛM , real numbers {xj }M
j=1 and some constant a > 0, we have


M
M
X
X
λj xj − aK(λ, π) ≤ a log 
πj exj /a  ,
j=1

j=1

where the equality is obtained when (xj /a) − log(λj /πj ) is a constant for 1 ≤ j ≤ M .
Proof: The result follows directly from Jensen’s Inequality as


M
M
M
X
X
X
exp 
λj ((xj /a) − log(λj /πj )) ≤
λj exp ((xj /a) − log(λj /πj )) =
πj exj /a .
j=1

j=1

j=1

Now by setting xj = −νkfj − hk22 and a = 2ω 2 in Lemma 3, we obtain


M
X
ν
min ν
λj kfj − hk22 + 2ω 2 K(λ, π) −
kh − Y k22
M
1−ν
λ∈Λ
j=1


M
X
ν
2
2
=−
kh − Y k22 − 2ω 2 log 
πj e−νkfj −hk2 /2ω  ,
1−ν
j=1

which implies that
T (h) = min S(λ, h).
λ∈ΛM

21

In addition, it is easy to verify that
Q(λ) = maxn S(λ, h),
h∈R

where the minimum is achieved at h = fλ .
Now let ĥ be the maximizer of T (h) in (21), then by setting the derivative of (21) to zero, it
is easy to observe that there exists a corresponding λ̂ so that (λ̂, ĥ) ∈ A ∩ B. This means that
A ∩ B 6= ∅.
Now consider any (λ0 , h0 ) ∈ A ∩ B. We have
Q(λ0 ) ≥ min Q(λ) = min maxn S(λ, h) ≥ maxn min S(λ, h).
λ∈ΛM

λ∈ΛM h∈R

h∈R λ∈ΛM

The third inequality is the well-known weak duality (e.g., Lemma 36.1 in [16]).
Also we have
max min S(λ, h) = maxn T (h) = T (ĥ) ≥ T (h0 ) .

h∈Rn λ∈ΛM

h∈R

We thus have
Q(λ0 ) ≥ min Q(λ) = min maxn S(λ, h) ≥ maxn min S(λ, h) = maxn T (h) ≥ T (h0 ).
λ∈ΛM

λ∈ΛM h∈R

h∈R λ∈ΛM

h∈R

Our target is now to prove Q(λ0 ) = T (h0 ). Since (λ0 , h0 ) ∈ A ∩ B, we have

1
1−ν
0


 h = ν Y − ν fλ0 ,

exp − 2ων 2 kfj − h0 k22 πj

0

 .
 λj = PM
0 2
ν
i=1 exp − 2ω 2 kfi − h k2 πi
It follows that for all j:
M
X
i=1



exp − 2ων 2 kfj − h0 k22 πj
ν
0 2
exp − 2 kfi − h k2 πi =
,
2ω
λ0j


which implies that
log

M
X
i=1


ν
exp − 2 kfi − h0 k22 πi
2ω


!
=−

=

ν
kfj − h0 k22 − log(λ0j /πj )
2ω 2

M
X
i=1

 ν

λ0i − 2 kfi − h0 k22 − log(λ0i /πi ) ,
2ω

where the second equation is from summing up two sides of the first equation with weight λ0i over
i = 1, . . . , M .

22

Plug back into T (h0 ), we obtain
#
" M

X 
ν
ν
T (h0 ) = −
kh0 − Y k22 − 2ω 2
λ0i − 2 kfi − h0 k22 − log(λ0i /πi )
1−ν
2ω
i=1

=−

ν
kh0 − Y k22 + ν
1−ν

= kfλ0 − Y k22 + ν

M
X

M
X

λ0i kfi − h0 k22 + 2ω 2 K(λ0 , π)

i=1

λ0i kfi − fλ0 k22 + 2ω 2 K(λ0 , π)

i=1
0

= Q(λ ) ,
where the third equality is obtained by plugging in h0 = ν1 Y −

1−ν
ν fλ0 .

Therefore,

Q(λ0 ) = min Q(λ) = min maxn S(λ, h) = maxn min S(λ, h) = maxn T (h) = T (h0 ).
λ∈ΛM

λ∈ΛM h∈R

h∈R λ∈ΛM

h∈R

Since Q(·) is strictly convex and T (·) is strictly concave, we have h0 = ĥ is the unique solution
of maxh T (h), and λ0 = λ is the unique solution of minλ Q(λ). Using h0 = ν1 Y − 1−ν
ν fλ0 , we have
ĥ =

1
1−ν
Y −
f Q.
ν
ν λ

This proves that A ∩ B contains the unique point (λQ , ĥ).

A.3

Proof of Proposition 1

The strong convexity of log J(·) in (14) implies that

2 
kψ̂ − ψ X (ω 2 , ν)k22 ≤
log J(ψ̂) − log J(ψ X (ω 2 , ν)) ≤ 2/A1 .
A1
Now plug the above inequality into the following equation
kψ̂ − ηk22 = kψ X (ω 2 , ν) − ηk22 + 2kψ̂ − ψ X (ω 2 , ν)k2 kψ X (ω 2 , ν) − ηk2 + kψ̂ − ψ X (ω 2 , ν)k22 ,
we obtain the desired bound.

A.4

Proof of Proposition 2

From definition, ψ X (ω 2 , ν) = fλ with λ ∈ ΛM defined as


1
1−ν
2
2
2
λj ∝ πj exp − 2 kfj − Y k2 +
kψ X (ω , ν) − fj k2 .
2ω
2ω 2
For any j = 1, . . . , M ,


log J(ψ (k) ) = log J ψ (k−1) + αk (fJ (k) − ψ (k−1) )


≤ log J ψ (k−1) + αk (fj − ψ (k−1) )
≤ log J(ψ (k−1) ) + αk (fj − ψ (k−1) )>

23

∇J(ψ (k−1) )
J(ψ (k−1) )

+ 2αk2 A3 ,

where the first inequality comes from definition, the second inequality is from Taylor expansion at
ψ (k−1) and (15) in Lemma 1 with the fact that kfj − ψ (k−1) k22 ≤ 4L2 .
We multiply the above inequality by λj and sum over j to obtain
log J(ψ

(k)

) ≤ log J(ψ

(k−1)

) + αk

M
X

λj (fj − ψ (k−1) )>

j=1

∇J(ψ (k−1) )
J(ψ

= log J(ψ (k−1) ) + αk (ψ X (ω 2 , ν) − ψ (k−1) )>

(k−1)

)

+ 2αk2 A3

∇J(ψ (k−1) )
J(ψ (k−1) )

+ 2αk2 A3

≤ log J(ψ (k−1) ) + αk (log J(ψ X (ω 2 , ν)) − log J(ψ (k−1) )) + 2αk2 A3 ,
where the last inequality follows from the convexity of log J(ψ).
Denote by δk = log J(ψ (k) ) − log J(ψ X (ω 2 , ν)), it follows that
δk ≤ (1 − αk )δk−1 + 2αk2 A3 .

P
We now bound δ0 . If we let µj ∝ πj exp − 2ω1 2 kfj − Y k22 such that M
j=1 µj = 1, then
δ0 = log J(ψ (0) ) − log J(ψ X (ω 2 , ν))




X
X
1 − ν (0)
1−ν
2
2
2
= log
µj exp
kψ − fj k2 − log
µj exp
kψ X (ω , ν) − fj k2
2ω 2
2ω 2
j
j




M
X
1
−
ν
≤ log 
µj exp
kψ (0) − fj k22 
2ω 2
j=1

1−ν 2
≤
L ≤ 2A3 .
2ω 2

(48)

The claim thus hold for δ0 . By mathematical induction, if δk−1 ≤
δk ≤ (1 − αk )δk−1 + 2αk2 A3 ≤ (1 − 2/(k + 1))

8A3
k+2 ,

then

8A3
8A3
+ 2(2/(k + 1))2 A3 ≤
.
k+2
k+3

This proves the desired bound.

A.5

Proof of Proposition 3

We have
kψ (k) − ψ X (ω 2 , ν)k22 ≤


2 8A3
16A3
2 
log J(ψ (k) ) − log J(ψ X (ω 2 , ν)) ≤
=
,
A1
A1 k + 3
A1 (k + 3)

where the first inequality comes from Taylor expansion at point ψ X (ω 2 , ν), with using (14) in
Lemma 1 and ∇J(ψ 2 ); and the second inequality is from Proposition 2. It follows that
kψ (k) − ηk22 = k(ψ (k) − ψ X (ω 2 , ν)) + (ψ X (ω 2 , ν) − η)k22
≤ kψ X (ω 2 , ν) − ηk22 + 2kψ X (ω 2 , ν) − ηk2 kψ (k) − ψ X (ω 2 , ν)k2 + kψ (k) − ψ X (ω 2 , ν)k22
s
16A3
16A3
≤ kψ X (ω 2 , ν) − ηk22 + 2
kψ X (ω 2 , ν) − ηk2 +
.
A1 (k + 3)
A1 (k + 3)

24

A.6

Proof of Proposition 4

The update operation implies
ψ (k) = ψ (k−1) − tk ∇ log J(ψ (k−1) ).
Then,
log J(ψ (k) ) = log J(ψ (k−1) − tk ∇ log J(ψ (k−1) ))
≤ log J(ψ (k−1) ) − tk k∇ log J(ψ (k−1) ))k22 + (A2 /2)t2k k∇ log J(ψ (k−1) ))k22
= log J(ψ (k−1) ) − (tk − (A2 /2)t2k )k∇ log J(ψ (k−1) ))k22 ,
where the inequality is from (15). By subtracting log J(ψ X (ω 2 , ν)) by each side, we have
log J(ψ (k) ) − log J(ψ X (ω 2 , ν)) ≤ log J(ψ (k−1) ) − log J(ψ X (ω 2 , ν))
− (tk − (A2 /2)t2k )k∇ log J(ψ (k−1) ))k22 .

(49)

Also from (15) we have


k∇ log J(ψ (k−1) ))k22 ≥ 2A1 log J(ψ (k−1) ) − log J(ψ X (ω 2 , ν)) .

(50)

Choose fixed step size tk = s ∈ (0, 2/A2 ) for any k > 0. Combining (49) and (50) results


log J(ψ (k) ) − log J(ψ X (ω 2 , ν)) ≤ [1 − 2A1 (s − (A2 /2)s2 )] log J(ψ (k−1) ) − log J(ψ X (ω 2 , ν)) .
It follows that


log J(ψ (k) ) − log J(ψ X (ω 2 , ν)) ≤ [1 − 2A1 (s − (A2 /2)s2 )]k log J(ψ (0) ) − log J(ψ X (ω 2 , ν)) .

A.7

Proof of Proposition 5

We choose tk = s as in Remark 1. Then,

2 
log J(ψ (k) ) − log J(ψ X (ω 2 , ν))
A1
2
≤
(1 − A1 /A2 )k log J(ψ (0) )
A1
2
1−ν 2
(1 − A1 /A2 )k
≤
L
A1
2ω 2
= L2 (1 − A1 /A2 )k ,

kψ (k) − ψ X (ω 2 , ν)k22 ≤

where the first inequality comes from Taylor expansion at point ψ X (ω 2 , ν), with using (14) in
Lemma 1 and ∇ log J(ψ X (ω 2 , ν)) = 0; the second inequality is from Proposition 4; and the third
inequality is from assumption (11) resulting log J(ψ (0) ) ≤ 1−ν
L2 . It follows that
2ω 2
kψ (k) − ηk22 = k(ψ (k) − ψ X (ω 2 , ν)) + (ψ X (ω 2 , ν) − η)k22
≤ kψ X (ω 2 , ν) − ηk22 + 2kψ X (ω 2 , ν) − ηk2 kψ (k) − ψ X (ω 2 , ν)k2 + kψ (k) − ψ X (ω 2 , ν)k22
q
≤ kψ X (ω 2 , ν) − ηk22 + 2 L2 (1 − A1 /A2 )k kψ X (ω 2 , ν) − ηk2 + L2 (1 − A1 /A2 )k .

25

A.8

Proof of Proposition 6

Given Y , the expectation is with respect to the randomness from the MH algorithm. For k > 0,
P
(k−1)
(k−1)
fj . Then in Algorithm 2, we update
uT
from Algorithm 3 is an estimator of fλ(k−1) = M
j=1 λj
(k)
ψ by
1−ν
(k−1)
ψ (k) = ψ (k−1) − tk 2 (ψ (k−1) − uT
).
ω
Denote v (k−1) =

1−ν
(ψ (k−1)
ω2

(k−1)

− uT

E[v (k−1) |ψ (k−1) ] =

), then we have

1 − ν (k−1)
(ψ
− fλ(k−1) ) = ∇ log J(ψ (k−1) ),
ω2

and
kCOV [v

(k−1)

|ψ

(k−1)


]kop =

1−ν
ω2

2
kCOV

(k−1)
[uT
|ψ (k−1) ]kop


≤

1−ν
ω2

2

s2 .

It follows that
log J(ψ (k) ) = log J(ψ (k−1) − tk v (k−1) )
≤ log J(ψ (k−1) ) − tk ∇ log J(ψ (k−1) )> v (k−1) + (A2 /2)t2k kv (k−1) k22 ,
where the inequality is from (15). Then, by subtracting log J(ψ X (ω 2 , ν)) from each side of the
above equation and taking expectation conditioned on ψ (k−1) , we have
E[δk |ψ (k−1) ] ≤ δk−1 − tk ∇ log J(ψ (k−1) )> E[v (k−1) |ψ (k−1) ] + (A2 /2)t2k E[kv (k−1) k22 |ψ (k−1) ]
!


1−ν 2 2
(k−1) 2
(k−1) 2
2
)k2 + n
≤ δk−1 − tk k∇ log J(ψ
)k2 + (A2 /2)tk k∇ log J(ψ
s
ω2


1−ν 2 2
1
1
(k−1) 2
= δk−1 −
ns ,
k∇ log J(ψ
)k2 +
2A2
2A2
ω2
where δk = log J(ψ (k) ) − log J(ψ X (ω 2 , ν)) and tk = s = 1/A2 as in Remark 1. Combining the
above inequality with


(51)
k∇ log J(ψ (k−1) ))k22 ≥ 2A1 log J(ψ (k−1) ) − log J(ψ X (ω 2 , ν)) ,
which is from (15), we have
E[δk |ψ (k−1) ] ≤ δk−1 (1 − A1 /A2 ) +

A21 2
ns .
2A2

Therefore, it follows that
E[δk ] ≤ E[δk−1 ](1 − A1 /A2 ) +

A21 2
ns ,
2A2

E[δk ] ≤ E[δ0 ](1 − A1 /A2 )k +

A1 2
ns .
2

and we have

26

References
[1] Dong Dai, Philippe Rigollet, Tong Zhang, et al. Deviation optimal learning using greedy
q-aggregation. Ann. Statist., 40(3):1878–1905, 2012.
[2] Alexandre B Tsybakov. Optimal rates of aggregation. In Learning Theory and Kernel Machines
(Lecture Notes in Computer Science), volume 2777, pages 303–313. Springer, 2003.
[3] Philippe Rigollet et al. Kullback–leibler aggregation and misspecified generalized linear models.
Ann. Statist., 40(2):639–665, 2012.
[4] Philippe Rigollet, Alexandre B Tsybakov, et al. Sparse estimation by exponential weighting.
Statistical Science, 27(4):558–575, 2012.
[5] Jean-Yves Audibert. Progressive mixture rules are deviation suboptimal. In Advances in
Neural Information Processing Systems, pages 41–48. MIT Press, 2008.
[6] Guillaume Lecué and Shahar Mendelson. Aggregation via empirical risk minimization. Probab.
Theory Rel. Fields, 145(3):591–613, 2009.
[7] Stéphane Gaı̈ffas and Guillaume Lecué. Hyper-sparse optimal aggregation. J. Mach. Learn.
Res., 12(Jun):1813–1833, 2011.
[8] Dong Dai and Tong Zhang. Greedy model averaging. In Advances in Neural Information
Processing Systems, pages 1242–1250. MIT Press, 2011.
[9] Arnak Dalalyan and Alexandre Tsybakov. Aggregation by exponential weighting and sharp
oracle inequalities. Learning Theory (Lecture Notes in Computer Science), pages 97–111, 2007.
[10] Arnak Dalalyan and Alexandre B Tsybakov. Aggregation by exponential weighting, sharp
pac-bayesian bounds and sparsity. Mach. Learn., 72(1):39–61, 2008.
[11] Gilbert Leung and Andrew R Barron. Information theory and mixing least-squares regressions.
IEEE Trans. Inf. Theory, 52(8):3396–3410, 2006.
[12] Marguerite Frank and Philip Wolfe. An algorithm for quadratic programming. Naval Res.
Logistics Quart., 3(1-2):95–110, 1956.
[13] Lee K. Jones. A simple lemma on greedy approximation in Hilbert space and convergence rates
for projection pursuit regression and neural network training. Ann. Statist., 20(1):608–613,
1992.
[14] Andrew R. Barron. Universal approximation bounds for superposition of a sigmoidal function.
IEEE Trans. Inf. Theory, 39(3):930–945, 1993.
[15] Jean-Michel Marin and Christian Robert. Bayesian core: a practical approach to computational
Bayesian statistics. Springer, 2007.
[16] Ralph Tyrell Rockafellar. Convex analysis. Princeton Univ. Press, 2015.

27

