DISTRIBUTIONALLY ROBUST GROUPWISE REGULARIZATION
ESTIMATOR

arXiv:1705.04241v1 [] 11 May 2017

BLANCHET, J. AND KANG, Y.

Abstract. Regularized estimators in the context of group variables have been applied
successfully in model and feature selection in order to preserve interpretability. We formulate
a Distributionally Robust Optimization (DRO) problem which recovers popular estimators,
such as Group Square Root Lasso (GSRL). Our DRO formulation allows us to interpret
GSRL as a game, in which we learn a regression parameter while an adversary chooses a
perturbation of the data. We wish to pick the parameter to minimize the expected loss under
any plausible model chosen by the adversary - who, on the other hand, wishes to increase
the expected loss. The regularization parameter turns out to be precisely determined by the
amount of perturbation on the training data allowed by the adversary. In this paper, we
introduce a data-driven (statistical) criterion for the optimal choice of regularization, which
we evaluate asymptotically, in closed form, as the size of the training set increases. Our
easy-to-evaluate regularization formula is compared against cross-validation, showing good
(sometimes superior) performance.

1. Introduction
Group Lasso (GR-Lasso) estimator is a generalization of the Lasso estimator (see Tibshirani (1996)). The method focuses on variable selection in settings where some predictive
variables, if selected, must be chosen as a group. For example, in the context of the use of
dummy variables to encode a categorical predictor, the application of the standard Lasso
procedure might result in the algorithm including only a few of the variables but not all of
them, which could make the resulting model difficult to interpret. Another example, where
the GR-Lasso estimator is particularly useful, arises in the context of feature selection. Once
again, a particular feature might be represented by several variables, which often should be
considered as a group in the variable selection process.
The GR-Lasso estimator was initially developed for the linear regression case (see Yuan and
Lin (2006)), but a similar group-wise regularization was also applied to logistic regression in
Meier et al. (2008). A brief summary of GR-Lasso technique type of methods can be found
in Friedman et al. (2010).
Recently, Bunea et al. (2014) developed a variation of the GR-Lasso estimator, called the
Group-Square-Root-Lasso (GSRL) estimator, which is very similar to the GR-Lasso estimator. The GSRL is to the GR-Lasso estimator what sqrt-Lasso, introduced in Belloni et al.
(2011), is to the standard Lasso estimator. In particular, GSRL has a superior advantage
over GR-Lasso, namely, that the regularization parameter can be chosen independently from
Columbia University. New York, NY 10027, United States.
E-mail addresses: jose.blanchet@columbia.edu, yang.kang@columbia.edu.
Date: May 12, 2017.
Key words and phrases. Distributionally Robust Optimization, Group Lasso, Optimal Transport.
1

2

the standard deviation of the regression error in order to guarantee the statistical consistency
of the regression estimator (see Belloni et al. (2011), and Bunea et al. (2014)).
Our contribution in this paper is to provide a DRO representation for the GSRL estimator, which is rich in interpretability and which provides insights to optimally select (using
a natural criterion) the regularization parameter without the need of time-consuming crossvalidation. We compute the optimal regularization choice (based on a simple formula we
derive in this paper) and evaluate its performance empirically. We will show that our method
for the regularization parameter is comparable, and sometimes superior, to cross-validation.

In order to describe our contributions more precisely, let us briefly describe the GSRL estimator. We choose the context of linear regression to simplify the exposition, but an entirely
analogous discussion applies to the context of logistic regression.
Consider a given a set of training data {(X1 , Y1 ), . . . , (Xn , Yn )}. The input Xi ∈ Rd is a
vector of d predicting variables, and Yi ∈ R is the response variable. (Throughout the paper
any vector is understood to be a column vector and the transpose of x is denoted by xT .)
We use (X, Y ) to denote a generic sample from the training data set. It is postulated that
Yi = XiT β ∗ + ei ,
for some β ∗ ∈ Rd and errors {e1 , ..., en }. Under suitable statistical assumptions (such as
independence of the samples in the training data), one may be interested in estimating β ∗ .
2
Underlying, we consider the square loss function, i.e. l (x, y; β) = y − β T x , for the purpose
of this discussion but this choice, as we shall see, is not necessary.
Throughout the paper we will assume the following group structure for the space of predictors.
There are d¯ ≤ d mutually exclusive groups, which form a partition. More precisely, suppose
that G1 , . . . , Gd¯ satisfies that Gi ∩ Gj = ∅ for i 6= j, that G1 ∪ ... ∪ Gd¯ = {1, ..., d}, and the
Gi ’s are non-empty. We will use gi to denote the cardinality of Gi and shall write G for a
generic set in the partition and let g denote the cardinality of G.
We shall denote by x (G) ∈ Rg the sub-vector x ∈ Rd corresponding to G. So, if G =
T
{i1 , ..., ig }, then x (G) = Xi1 , . . . , Xig .
¯
¯ we define for each x ∈ Rd ,
Next, given p, s ≥ 1, and α ∈ Rd++ (i.e. αi > 0 for 1 ≤ i ≤ d)
!1/s
d¯
X
s
s
(1)
kxkα-(p,s) =
αi kx (Gi )kp
,
i=1

where kx (Gi )kp denotes the p-norm of x (Gi ) in Rgi . (We will study fundamental properties
of kxkα-(p,s) as a norm in Proposition 1.)
Let Pn be the empirical distribution function, namely,
n

Pn (dx, dy) :=

1X
δ{(Xi ,Yi )} (dx, dy).
n i=1

Throughout out the paper we use the notation EP [·] to denote expectation with respect to a
probability distribution P .

3

The GSRL estimator takes the form
v
u n


u1 X
1/2
min t
l (Xi, Yi ; β) + λ kβkg̃−1 −(2,1) = min EPn [l (X, Y ; β)] + λ kβk√g̃−(2,1) ,
β
β
n i=1
where λ is the so-called regularization parameter. The previous optimization problem can
be easily solved using standard convex optimization techniques as explained in Belloni et al.
(2011) and Bunea et al. (2014).
Our contributions in this paper can now be explicitly stated. We introduce a notion of
discrepancy, Dc (P, Pn ), discussed in Section 2, between Pn and any other probability measure
P , such that


1/2
1/2
(2)
min
max
EP [l (X, Y ; β)] = min EPn [l (X, Y ; β)] + δ 1/2 kβkα−(p,s) .
β

P :Dc (P,Pn )≤δ

β

Using this representation, which we formulate, together with its logistic regression analogue,
in Section 2.2.1 and Section 2.2.2, we are able to draw the following insights:
I) GSRL can be interpreted as a game in which we choose a parameter (i.e. β) and an
adversary chooses a “plausible” perturbation of the data (i.e. P ); the parameter δ controls
the degree in which Pn is allowed to be perturbed to produce P . The value of the game is
dictated by the expected loss, under EP , of the decision variable β.
II) The set Uδ (Pn ) = {P : Dc (P, Pn ) ≤ δ} denotes the set of distributional uncertainty. It
represents the set of plausible variations of the underlying probabilistic model which are
reasonably consistent with the data.
III) The DRO representation (2) exposes the role of the regularization parameter. In particular, because λ = δ 1/2 , we conclude that λ directly controls the size of the distributionally
uncertainty and should be interpreted as the parameter which dictates the degree to which
perturbations or variations of the available data should be considered.
IV) As a consequence of I) to III), the DRO representation (2) endows the GSRL estimator
with desirable generalization properties. The GSRL aims at choosing a parameter, β, which
should perform well for all possible probabilistic descriptions which are plausible given the
data.
Naturally, it is important to understand what types of variations or perturbations are measured by the discrepancy Dc (P, Pn ). For example, a popular notion of the discrepancy is
the Kullback-Leibler (KL) divergence. However, KL divergence has the limitation that only
considers probability distributions which are supported precisely on the available training
data, and therefore potentially ignores plausible variations of the data which could have an
adverse impact on generalization risk.
In the rest of the paper we answer the following questions. First, in Section 2 we explain the
nature of the discrepancy Dc (P, Pn ), which we choose as an Optimal Transport discrepancy.
We will see that Dc (P, Pn ) can be computed using a linear program.
Intuitively, Dc (P, Pn ) represents the minimal transportation cost for moving the mass encoded by Pn into a sinkhole which is represented by P . The cost of moving mass from location
u = (x, y) to w = (x0 , y 0 ) is encoded by a cost function c (u, w) which we shall discuss and
this will depend on the α-(p, s) norm that we defined in (1). The subindex c in Dc (P, Pn )

4

represents the dependence on the chosen cost function.
The next item of interest is the choice of δ, again the discussion of items I) to III) of the
DRO formulation (2) provides a natural way to optimally choose δ. The idea is that every
model P ∈ Uδ (Pn ) should intuitively represent a plausible variation of Pn and therefore
β P = arg min {EP [l (X, Y ; β)] : β} is a plausible estimate of β ∗ . The set {β P : P ∈ Uδ (Pn )}
therefore yields a confidence region for β ∗ which is increasing in size as δ increases. Hence,
it is natural to minimize δ to guarantee a target confidence level (say 95%). In Section 3 we
explain how this optimal choice can be asymptotically computed as n → ∞.
Finally, it is of interest to investigate if the optimal choice of δ (and thus of λ) actually
performs well in practice. We compare performance of our (asymptotically) optimal choice
of λ against cross-validation empirically in Section 4. We conclude that our choice is quite
comparable to cross validation.
Before we continue with the program that we have outlined, we wish to conclude this Introduction with a brief discussion of work related to the methods discussed in this paper.
Connections between regularized estimators and robust optimization formulations have been
studied in the literature. For example, the work of Xu et al. (2009) investigates determinist perturbations on the predictor variables to quantify uncertainty. In contrast, our DRO
approach quantifies perturbations from the empirical measure. This distinction allows us to
statistical theory which is key to optimize the size of the uncertainty, δ (and thus the regularization parameter) in a data-driven way. The work of Shafieezadeh-Abadeh et al. (2015)
provides connections to regularization in the setting of logistic regression in an approximate
form. More importantly, Shafieezadeh-Abadeh et al. (2015) propose a data-driven way to
choose the size of uncertainty, δ, which is based on the concentration of measure results. The
concentration of measure method depends on restrictive assumptions and leads to suboptimal choices, which deteriorate poorly in high dimensions.
The present work is a continuation of the line of research development in Blanchet et al.
(2016), which concentrates only on classical regularized estimators without the group structure. Our current contributions require the development of duality results behind the α-(p, t)
norm which closely parallels that of the standard duality between lp and lq spaces (with
1/p + 1/q = 1) and a number of adaptations and interpretations that are special to the group
setting only.
2. Optimal Transport and DRO
2.1. Defining the optimal transport discrepancy. Let c : Rd+1 × Rd+1 → [0, ∞] be
lower semicontinuous and we assume that c(u, w) = 0 if and only if u = w. For reasons that
will become apparent in the sequel, we will refer to c (·) as a cost function.
Given two distributions P and Q, with supports SP ⊆ Rd+1 and SQ ⊆ Rd+1 , respectively, we
define the optimal transport discrepancy, Dc , via
(3)

Dc (P, Q) = inf {Eπ [c (U, W )] : π ∈ P (SP × SQ ) , πU = P, πW = Q},
π

where P (SP × SQ ) is the set of probability distributions π supported on SP × SQ , and πU
and πW denote the marginals of U and W under π, respectively.
If, in addition, c (·) is symmetric (i.e. c (u, w) = c (w, u)), and there exists % ≥ 1 such that
c1/% (u, w) ≤ c1/% (u, v) + c1/% (v, w) (i.e. c1/% (·) satisfies the triangle inequality), it can be

5
1/%

easily verified (see Villani (2008)) that Dc (P, Q) is a metric on the probability measures.
For example, if c (u, w) = ku − wk%q for q ≥ 1 (where ku − wkq denotes the lq norm in Rd+1 )
then Dc (·) is known as the Wasserstein distance of order %.
Observe that (3) is obtained by solving a linear programming problem. For example, suppose
that Q = Pn , so SPn = {(Xi , Yi )}ni=1 , and let P be supported in some finite set SP then,
using U = (X, Y ), we have that Dc (P, Pn ) is obtained by computing
X X
(4)
min
c (u, w) π (u, w)
π

s.t.

u∈SP w∈SPn

X

π (u, w) =

u∈SP

X

1
∀ w ∈ SPn
n

π (u, w) = P ({u}) ∀ u ∈ SP

w∈SPn

π (u, w) ≥ 0 ∀ (u, w) ∈ SP × SPn .
A completely analogous linear program (LP), albeit an infinite dimensional one, can be
defined if SP has infinitely many elements. This LP has been extensively studied in great
generality in the context of Optimal Transport under the name of Kantorovich’s problem
(see Villani (2008))).
Note that Kantorovich’s problem is always feasible (take π with independent marginals, for
example). Moreover, under our assumptions on c, if the optimal value is finite, then there is
an optimal solution π ∗ (see Chapter 1 of Villani (2008))).
It is clear from the formulation of the previous LP that Dc (P, Pn ) can be interpreted as
the minimal cost of transporting mass from Pn to P , assuming that the marginal cost of
transporting the mass from u ∈ SP to w ∈ SPn is c (u, w). It is also not difficult to realize
from the assumption that c (u, w) = 0 if and only if u = w that Dc (P, Pn ) = 0 if and
only if P = Pn . We shall discuss, for instance, how to choose c (·) to recover (2) and the
corresponding logistic regression formulation of GR-Lasso.
2.2. DRO Representation of GSRL Estimators. In this section, we will construct a
cost function c (·) to obtain the GSRL (or GR-Lasso) estimators. We will follow an approach
introduced in Blanchet et al. (2016) for the context of square-root Lasso (SR-Lasso) and
regularized logistic regression estimators.
2.2.1. GSRL Estimators for Linear Regression. We start by assuming precisely the linear
regression setup described in the Introduction and leading to (2). Given α = (α1 , ..., αd¯)T ∈
T
d¯
R++
define α−1 = α1−1 , ..., αd−1
. Now, underlying there is a partition G1 , ..., Gd¯ of {1, ..., d}
¯
and given q, t ∈ [1, ∞] we introduce the cost function

kx − x0 k%α−1 -(q,t) if y = y 0
0 0
(5)
c ((x, y) , (x , y )) =
,
∞
if y 6= y 0
where, following (1), we have that
kx −

%
x0 kα−1 -(q,t)

=

d¯
X
i=1

!%/t
αi−t kx (Gi ) −

t
x0 (Gi )kq

.

6

Then, we obtain the following result.
Theorem 1 (DRO Representation for Linear Regression GSRL). Suppose that q, t ∈ [1, ∞]
d¯
are given and c (·) is defined as in (5) for % = 2. Then, if l (x, y; β) =
and α ∈ R++
2
T
y − x β we obtain
√
min
sup
(EP [l (X, Y ; β)])1/2 = min (EPn [l (X, Y ; β)])1/2 + δ kβkα-(p,s) ,
β∈Rd P :Dc (P,Pn )≤δ

β∈Rd

where 1/p + 1/q = 1, and 1/s + 1/t = 1.
√
¯ we end up
We remark that choosing p = q = 2, t = ∞, s = 1, and αi = gi for i ∈ {1, ..., d}
obtaining the GSRL estimator formulated in Bunea et al. (2014)).
We note that the cost function c (·) only allows mass transportation on the predictors (i.e
X), but no mass transportation is allowed on the response variable Y . This implies that
the GSRL estimator implicitly assumes that distributional uncertainty is only present on
prediction variables (i.e. variations on the data only occurs through the predictors).
2.2.2. GR-Lasso Estimators for Logistic Regression. We now discuss GR-Lasso for classification problems. We consider a training data set of the form {(X1 , Y1 ), . . . , (Xn , Yn )}. Once
again, the input Xi ∈ Rd is a vector of d predictor variables, but now the response variable
Yi ∈ {−1, 1} is a categorical variable. In this section we shall consider as our loss function
the log-exponential function, namely,

(6)
l (x, y; β) = log 1 + exp −yβ T x .
This loss function is motivated by a logistic regression model which we shall review in the sequel. But for the DRO representation formulation it is not necessary to impose any statistical
assumption. We then obtain the following theorem.
Theorem 2 (DRO Representation for Logistic Regression GR-Lasso). Suppose that q, t ∈
d¯
[1, ∞] and α ∈ R++
are given and c (·) is defined as in (5) for % = 1. Then, if l (x, y; β) is
defined as in (6) we obtain
min

sup

β∈Rd P :Dc (P,Pn )≤δ

EP [l (X, Y ; β)] = min EPn (l (X, Y ; β)) + δ kβkα-(p,s) ,
β∈Rd

where 1 ≤ q, t ≤ ∞, 1/p + 1/q = 1 and 1/s + 1/t = 1.
√
¯ and λ = δ we
We note that by taking p = q = 2, t = ∞, s = 1, αi = gi for i ∈ {1, ..., d},
recover the GR-Lasso logistic regression estimator from Meier et al. (2008).
As discussed in the previous subsection, the choice of c (·) implies that the GR-Lasso estimator implicitly assumes that distributionally uncertainty is only present on prediction
variables.
3. Optimal Choice of Regularization Parameter
Let us now discuss the mathematical formulation of the optimal criterion that we discussed
for choosing δ (and therefore the regularization parameter λ). We define
Λδ (Pn ) = {β P : P ∈ Uδ (Pn )},
as discussed in the Introduction, Λδ (Pn ) is a natural confidence region for β ∗ because each
element P in the distributional uncertainty set Uδ (Pn ) can be interpreted as a plausible

7

variation of the empirical data Pn . Then, given a confidence level 1 − χ (say 1 − χ = .95) we
wish to choose
δn∗ = inf{δ : P (β ∗ ∈ Λδ (Pn )) > 1 − χ}.
Note that in the evaluation of P (β ∗ ∈ Λδ (Pn )) the random element is Pn . So, we shall impose
natural probabilistic assumptions on the data generating process in order to asymptotically
evaluate δn∗ as n → ∞.
3.1. The Robust Wasserstein Profile Function. In order to asymptotically evaluate δn∗
we must recall basic properties of the so-called Robust Wassertein Profile function (RWP
function) introduced in Blanchet et al. (2016).
Suppose for each (x, y), the loss function l (x, y; ·) is convex and differentiable, then under
natural moment assumptions which guarantee that expectations are well defined, we have
that for
P ∈ Uδ (Pn ) = {P : Dc (P, Pn ) ≤ δ},
the parameter β P must satisfy
(7)



EP ∇β l X, Y ; β P = 0.

Now, for any given β, let us define
M (β) = {P : EP [∇β l (X, Y ; β)] = 0} ,
which is the set of probability measures P , under which β is the optimal risk minimization
parameter. We would like to choose δ as small as possible so that
(8)

Uδ (Pn ) ∩ M (β ∗ ) 6= ∅

with probability at least 1 − χ. But note that (8) holds if and only if there exists P such
that Dc (P, Pn ) ≤ δ and EP [∇β l (X, Y ; β ∗ )] = 0.
The RWP function is defined
(9)

Rn (β) = min{Dc (P, Pn ) : EP [∇β l (X, Y ; β)] = 0}.

In view of our discussion following (8), it is immediate that β ∗ ∈ Λδ (Pn ) if and only if
Rn (β ∗ ) ≤ δ, which then implies that
δn∗ = inf{δ : P (Rn (β ∗ ) ≤ δ) > 1 − χ}.
Consequently, we conclude that δn∗ can be evaluated asymptotically in terms of the 1 − χ
quantile of Rn (β ∗ ) and therefore we must identify the asymptotic distribution of Rn (β ∗ ) as
n → ∞. We illustrate intuitively the role of the RWP function and M (β) in Figure 1, where
RWP function Rn (β ∗ ) could be interpreted as the discrepancy distance between empirical
measure Pn and the manifold M (β ∗ ) associated with β ∗ .
Typically, under assumptions supporting the underlying model (as in the generalized linear
setting we considered), we will have that β ∗ is characterized by the estimating equation
(7). Therefore, under natural statistical assumptions one should expect that Rn (β ∗ ) → 0 as
n → ∞ at a certain rate and therefore δn∗ → 0 at a certain (optimal) rate. This then yields
an optimal rate of convergence to zero for the underlying regularization parameter. The next
subsections will investigate the precise rate of convergence analysis of δn∗ .

8

Figure 1. Intuitive Plot for the RWP function Rn (β) and the set M (β).
3.2. Optimal Regularization for GSRL Linear Regression. We assume, for simplicity, that the training data set {(X1 , Y1 ), . . . , (Xn , Yn )} is i.i.d. and that the linear relationship Yi = β ∗ T Xi + ei , holds with the errors {e1 , ..., en } being i.i.d. and independent of
{X1 , . . . , Xn }. Moreover, we assume that both the entries of Xi and the errors have finite
second moment and the errors have zero mean. 
2
Since in our current setting l (x, y; β) = y − xT β , then the RWP function (9) for linear
regression model is given as,



(10)
Rn (β) = min Dc (P, Pn ) : EP X Y − X T β = 0 .
P

Theorem 3 (RWP Function Asymptotic Results: Linear Regression). Under the assumptions imposed in this subsection and the cost function as given in (5), with % = 2,
n
h
io
 ∗ 2
∗
T
T
nRn (β ) ⇒ L1 := max 2σζ Z − E eζ − ζ X β α-(p,s) ,
ζ∈Rd

as n → ∞, where ⇒ means convergence in distribution and Z ∼ N (0, Σ) with Σ = V ar(X).
Moreover, we can observe the more tractable stochastic upper bound,
D
E [e2 ]
kZk2α−1 -(q,t) .
L1 ≤ L2 :=
E [e2 ] − (E [|e|])2
We now explain how to use Theorem 3 to set the regularization parameter in
GSRL linear regression:
(1) Estimate the 1 − χ quantile of kZk2α−1 -(q,t) . We use use η̂1−χ to denote the estimator
for this quantile. This step involves estimating Σ from the training data.
(2) The regularization parameter λ in the GSRL linear regression takes the form
√
−1/2
1/2
λ = δ = η̂1−χ n(1 − (E |e|)2 /Ee2 )
.
Note that the denominator in the previous expression must be estimated from the
training data.
Note that the regularization parameter for GSRL for linear regression chosen via our RWPI
asymptotic result does not depends on the magnitude of error e (see also the discussion in
Bunea et al. (2014)).

9

It is also possible to formulate the optimal regularization results for high-dimension setting,
where the number of predictors growth with sample size. We discuss the results in the
Appendix namely Section B.3.
3.3. Optimal Regularization for GR-Lasso Logistic Regression. We assume that the
training data set {(X1 , Y1 ), . . . , (Xn , Yn )} is i.i.d.. In addition, we assume that the Xi ’s have
a finite second moment and also that they possess a density with respect to the Lebesgue
measure. Moreover, we assume a logistic regression model; namely,

(11)
P (Yi = 1|Xi ) = 1/ 1 + exp −XiT β ∗ ,
and P (Yi = −1|Xi ) = 1 − P (Yi = 1|Xi ).
In the logistic regression setting, we consider the log-exponential loss defined in (6). Therefore, the RWP function, (9), for logistic regression is




YX
(12)
Rn (β) = min Dc (P, Pn ) : EP
=0 .
1 + exp (Y X T β)
Theorem 4 (RWP Function Asymptotic Results: Logistic Regression). Under the assumptions imposed in this subsection and the cost function as given in (5) with % = 1,
√
nRn (β ∗ ) ⇒ L3 := sup ζ T Z,
ζ∈A

as n → ∞, where

Z∼N

XX T
0, E
(1 + exp (Y X T β ∗ ))2




and
A=




T

ζ ∈ Rd : ess sup ζ T



y 1 + exp Y X β

∗



Id×d − XX

(1 + exp (Y X T β ∗ ))2

X,Y

T

α-(p,s)



≤1 .


Further, the limit law L3 follows the simpler stochastic bound,
D

L3 ≤ L4 := Z̃

α−1 -(q,t)

,

where Z̃ ∼ N (0, Σ).
We now explain how to use Theorem 4 to set the regularization parameter in
GR-Lasso logistic regression.
(1) Estimate the 1 − χ quantile of L4 . We use use η̂1−χ to denote the estimator for this
quantile. This step involves estimating Σ from the training data.
(2) We choose the regularization parameter λ in the GR-Lasso problem to be,
√
λ = δ = η̂1−χ / n.

10

4. Numerical Experiments
We proceed to numerical experiments on both simulated and real data to verify the performance of our method for choosing the regularization parameter. We apply “grpreg” in R,
from Breheny and Breheny (2016), to solve GR-Lasso for logistic regression. For GSRL for
linear regression, we consider apply the “grpreg” solver for the GR-Lasso problem combined
with the iterative procedure discussed in Section 2 of Sun and Zhang (2011) (see also Section
5 of Li et al. (2015) for the Lasso counterpart of such numerical procedure).
Data preparation for simulated experiments: We borrow the setting from example III
in Yuan and Lin (2006), where the group structure is determined by the third order polynomial expansion. More specifically, we assume that we have 17 random variables Z1 , . . . , Z16
and W , they are i.i.d. and
√ follow the normal distribution. The covariates X1 , . . . , X16 are
given as Xi = (Zi + W ) / 2. For the predictors, we consider each covariate and its second
and third order polynomial, i.e. Xi , Xi2 and Xi3 . In total, we have 48 predictors.
For linear regression: The response Y is given by
Y = β3,1 X3 + β3,2 X32 + β3,3 X33 + β5,1 X5 + β5,2 X52 + β5,3 X53 + e,
where β(·,·) coefficients draw randomly and e represents an independent random error.
For classification: We consider Y simulated by a Bernoulli distribution, i.e.


Y ∼ Ber 1/ 1 + exp − β3,1 X3 + β3,2 X32 + β3,3 X33 + β5,1 X5 + β5,2 X52 + β5,3 X53
.
We compare the following methods for linear regression and logistic regression: 1) groupwise
regularization with asymptotic results (in Theorem 3, 4) selected tuning parameter (RWPI
GRSL and RWPI GR-Lasso), 2) groupwise regularization with cross-validation (CV GRSL
and CV GR-Lasso), and 3) ordinary least square and logistic regression (OLS and LR).
We report the error as the square loss for linear regression and log-exponential loss for logistic
regression. The training error is calculated via the training data. The size of the training
data is taken to be n = 50, 100, 500 and 1000. The testing error is evaluated using a simulated
data set of size 1000 using the same data generating process described earlier. The mean
and standard deviation of the error are reported via 200 independent runs of the whole
experiment, for each sample size n.
The detailed results are summarized in Table 1 for linear regression and Table 2 for logistic
regression. We can see that our procedure is very comparable to cross validation, but it is
significantly less time consuming and all of the data can be directly used to estimate the
model parameter, by-passing significant data usage in the estimation of the regularization
parameter via cross validation
RWPI
Sample Size Training
n = 50
5.64 ± 1.16
n = 100
4.67 ± 0.70
n = 500
4.09 ± 0.29
n = 1000
4.02 ± 0.19

GSRL
CV GSRL
Testing
Training
Testing
9.15 ± 3.58 3.18 ± 1.07 7.66 ± 2.69
5.83 ± 1.38 3.61 ± 0.74 5.22 ± 1.05
4.16 ± 0.27 3.93 ± 0.3 4.12 ± 0.27
4.11 ± 0.26 3.95 ± 0.19 4.11 ± 0.26

OLS
Training
Testing
0.07 ± 0.09 80.98 ± 30.53
2.09 ± 0.44 73.35 ± 16.51
3.63 ± 0.27 73.08 ± 10.40
3.82 ± 0.19 72.28 ± 8.05

Table 1. Linear Regression Simulation Results.

11

Sample Size
n = 50
n = 100
n = 500
n = 1000

RWPI GR-Lasso
Training
Testing
.683 ± .016 .702 ± .014
.593 ± .038 .618 ± .029
.513 ± .021 .518 ± .019
.492 ± .016 .488 ± .017

CV GR-Lasso
Training
Testing
.459 ± .118 .628 ± .099
.450 ± .061 .551 ± .037
.461 ± .025 .493 ± .018
.491 ± .017 .488 ± .019

Logistic
Training
.002 ± .001
.042 ± .041
.083 ± .057
.442 ± .018

Regression
Testing
5.288 ± 1.741
4.571 ± 1.546
1.553 ± .355
.510 ± .028

Table 2. Logistic Regression Simulation Results.

We also validated our method using the Breast Cancer classification problem with data from
the UCI machine learning database discussed in Lichman (2013). The data set contains 569
samples with one binary response and 30 predictors. We consider all the predictors and their
first, second, and third order polynomial expansion. Thus, we end up having 90 predictors
divided into 30 groups. For each iteration, we randomly split the data into a training set
with 112 samples and the rest in the testing set. We repeat the experiment 500 times to
observe the log-exponential loss function for the training and testing error. We compare our
asymptotic results based GR-Lasso logistic regression (RWPI GR-Lasso), cross-validation
based GR-Lasso logistic regression (CV GR-Lasso), vanilla logistic regression (LR), and
regularized logistic regression (LRL1). We can observe, even when the sample size is small
as in the example, our method still provides very comparable results (see in Table 3).
LR
Training
0.0 ± 0.0

Testing
15.267 ± 5.367

LRL1
Training
Testing
.510 ± .215 .414 ± .173

RWPI GR-Lasso
CV GR-Lasso
Training
Testing
Training
Testing
.186 ± .032 .240 ± .098 .198 ± .041 .213 ± .041

Table 3. Numerical results for breast cancer data set.

5. Conclusion and Extensions
Our discussion of GSRL as a DRO problem has exposed rich interpretations which we have
used to understand GSRL’s generalization properties by means of a game theoretic formulation. Moreover, our DRO representation also elucidates the crucial role of the regularization parameter in measuring the distributional uncertainty present in the data. Finally,
we obtained asymptotically valid formulas for optimal regularization parameters under a
criterion which is naturally motivated, once again, thanks to our DRO formulation. Our
easy-to-implement formulas are shown to perform well compared to (time-consuming) cross
validation.
We strongly believe that our discussion in this paper can be easily extended to a wide range
of machine learning estimators. We envision formulating the DRO problem considering different types of models and cost functions. We plan to investigate algorithms which solve
the DRO problem directly (even if no direct regularization representation, as the one we
considered here, exists). Moreover, it is natural to consider different types of cost functions
which might improve upon the simple choice which, as we have shown, implies the GSRL
estimator. Questions related to alternative choices of cost functions are also under current
research investigations, and our progress will be reported in the near future.

12

References
[1] Belloni, A., Chernozhukov, V., and Wang, L. (2011). Square-root Lasso: pivotal recovery
of sparse signals via conic programming. Biometrika, 98(4):791–806.
[2] Blanchet, J., Kang, Y., and Murthy, K. (2016). Robust Wasserstein profile inference and
applications to machine learning. arXiv preprint arXiv:1610.05627.
[3] Breheny, P. and Breheny, M. P. (2016). Package grpreg.
[4] Bunea, F., Lederer, J., and She, Y. (2014). The group square-root Lasso: Theoretical
properties and fast algorithms. IEEE Transactions on Information Theory, 60(2):1313–1325.
[5] Friedman, J., Hastie, T., and Tibshirani, R. (2010). A note on the group Lasso and a
sparse group Lasso. arXiv preprint arXiv:1001.0736.
[6] Li, X., Zhao, T., Yuan, X., and Liu, H. (2015). The Flare package for high dimensional
linear regression and precision matrix estimation in R. JMLR, 16:553–557.
[7] Lichman, M. (2013). UCI machine learning repository.
[8] Meier, L., Van De Geer, S., and Bühlmann, P. (2008). The group Lasso for logistic
regression. Journal of the Royal Statistical Society: Series B, 70(1):53–71.
[9] Shafieezadeh-Abadeh, S., Esfahani, P. M., and Kuhn, D. (2015). Distributionally robust
logistic regression. In Advances in NIPS 28, pages 1576–1584.
[10] Sun, T. and Zhang, C.-H. (2011). Scaled sparse linear regression. arXiv preprint
arXiv:1104.4595.
[11] Tibshirani, R. (1996). Regression shrinkage and selection via the Lasso. Journal of the
Royal Statistical Society. Series B, 58(1):267–288.
[12] Villani, C. (2008). Optimal Transport: old and new. Springer Science & Business Media.
[13] Xu, H., Caramanis, C., and Mannor, S. (2009). Robust regression and Lasso. In Advances
in NIPS 21, pages 1801–1808.
[14] Yuan, M. and Lin, Y. (2006). Model selection and estimation in regression with grouped
variables. Journal of the Royal Statistical Society: Series B, 68(1):49–67.

Appendix A. Technical Proofs
We will first derive some properties for α-(p, s) norm we defined in (1), then we move to the
proof for DRO problem in Section A.2 and the optimal selection of regularization parameter
in Section A.3. We will focus on the proof for linear regression and leave the part for logistic
regression, which follows the similar techniques, in the Appendix B.
A.1. Basic Properties of the α-(p, s) Norm . The following Proposition, which describes
basic properties of the α-(p, s) norm, will be very useful in our proofs.
Proposition 1. For α − (p, s) norm defined for Rd as in (1) and the notations therein, we
have the following properties:
I) The dual norm of α − (p, s) norm is α−1 -(q, t) norm, where α−1 = (1/α1 , . . . , 1/αd¯)T ,
1/p + 1/q = 1, and 1/s + 1/t = 1 (i.e. p, q are conjugate and s, t are conjugate).
II) The Hölder inequality holds for the α-(p, s) norm, i.e. for a, b ∈ Rd , we have,
aT b ≤ kakα-(p,s) kbkα−1 -(q,t) ,

13

where the equality holds if and only if sign(a(Gj )i ) = sign(b(Gj )i ) and
1
|αj a(Gj )i |
b(Gj )
αj

q/p−t/s
t/s
kbkα−1 -(q,t)
q

q/p

1
=
b(Gj )i
αj

.

is true for all j = 1, . . . , d¯ and i = 1, . . . , gj .
The triangle inequality holds, i.e. for a, b ∈ Rd and a 6= 0, we have
kakα-(p,s) + kbkα-(p,s) ≥ ka + bkα-(p,s) ,
where the equality holds if and only if, there exists nonnegative τ , such that τ a = b.
Proof of Proposition 1. We first proceed to prove II). Let us consider any a, b ∈ Rd . We can
assume a, b 6= 0, otherwise the claims are immediate. The inner product (or dot product) of
a and b an be written as:
" gj
#
" gj
#
d¯
d¯
X
X
X
X
|a(Gj )i | · |b(Gj )i | .
a(Gj )i b(Gj )i ≤
aT b =
j=1

j=1

i=1

i=1

The equality holds for the above inequality if and only if a(Gj )i and b(Gj )i shares the same
¯ we consider the term in the bracket,
sign. For each fixed j = 1, . . . , d,
gj
X

|a(Gj )i | · |b(Gj )i | =

gj
X

αj |a(Gj )i | · |b(Gj )i | /αj ≤ kαj a(Gj )kp ·

i=1

i=1

1
b(Gj )
αj .

.
q

The above inequality is due to Hölder’s inequality for p−norm and the equality holds if and
only if
q
q
1
p
p 1
b(Gj ) |αj a(Gj )i | = kαj a(Gj )kp
b(Gj )i ,
αj .
αj
q
¯ we have,
is true for all i = 1, gj . Combining the above result for each j = 1, . . . , d,
T

a b≤

d¯
X
j=1

kαj a(Gj )kp ·

1
b(Gj )
αj

≤ kakα-(p,s) · kbkα−1 -(q,t) ,
q

where the final inequality is due to Hölder inequality applied to the vectors

T

T
1
1
kbG1 kq , . . . ,
(13) ã = α1 ka(G1 )kp , . . . , αd¯ ka(Gd¯)kp , and b̃ =
kb(Gd¯)kq
.
α1
αd¯
This proves the Hölder type inequality stated in the theorem. We can further observe that
the final inequality becomes equality if and only if
kbktα−1 −(q,t) kαj a(Gj )ksp = kaksα−(p,s)

1
b(Gj )
αj

t

,
q

¯ Combining the conditions for equalities hold for each inequality,
holds for all j = 1, . . . , d.
we conclude condition II) in the statement of the proposition.
Now we proceed to prove I). Recall the definition of a dual norm, i.e.
kbk∗α-(p,s) =

sup
a:kakα-(p,s) =1

aT b

14

. Now, choose b ∈ Rd , b 6= 0, and let us take a satisfying, kakα−(p,s) = 1 and

a(Gj )i =

sign(b(Gj )i )
αj

1
b(Gj )i
αj
1
b(Gj )
αj

q/p−t/s
q

q/p

.
t/s
kbkα−1 -(q,t)

By part II), we have that
kbk∗α−(p,s) =

sup
a:kakα−(p,s) =1

aT b = kakα−(p,s) kbkα−1 −(q,t) = kbkα−1 −(q,t) .

Thus we proved part I). Finally, let us discuss the triangle inequality. For any a, b ∈ Rd and
a, b 6= 0 we have
kakα-(p,s) + kbkα-(p,s)
" d¯
#1/s " d¯
#1/s
X
X
=
αj ka(Gj )ksp
+
αj kb(Gj )ksp
j=1

"
≥

d¯
X

j=1



αj ka(Gj )ksp + kb(Gj )ksp



#1/s

j=1

"
≥

d¯
X

#1/s
αj ka(Gj ) +

b(Gj )ksp

j=1

= ka + bkα-(p,s) .
For the above derivation, the first equality is due to definition in (1), Second equality is
applying the triangle inequality of s-norm for ã and b̃ defined in (13), where the equality
holds if and only if, there exist positive number τ̃ , such that τ̃ ã = b̃. Third inequality
¯ where the
is due to triangle equality of p-norm to a(Gj ) and b(Gj ) for each j = 1, . . . , d,
equality holds if and only if, there exists nonnegative numbers τj , such that τj a(Gj ) = b(Gj ).
Combining the equality condition for second and third estimate above, we can conclude the
equality condition for the triangle inequality for α-(p, s) norm is if and only if there exists a
non-negative number τ , such that τ a = b.

A.2. Proof of DRO for Linear Regression .
Proof of Theorem 1. Let us apply the strong duality results, as in the Appendix of Blanchet
et al. (2016), for worst-case expected loss function, which is a semi-infinity linear programming problem, and write the worst-case loss as,
(
)
n
n
o
h
i
X


1
2
2
sup yi − uT β − γ kxi − uk2α−1 -(q,t)
.
sup
EP Y − X T β
= min γδ −
γ≥0
n
u
P :Dc (P,Pn )≤δ
i=1

15

For each i, let us consider the inner optimization problem over u. We can denote ∆ = u − xi
and ei = yi −xTi β for notation simplicity, then the i−th inner optimization problem becomes,
o
n
2
2
T
T
2
ei + sup ∆ β − 2ei ∆ β − γ k∆kα−1 −(q,t)
∆


!2
 X

X
2
2
=ei + sup
|∆j | |βj | + 2 |ei |
|∆j | |βj | − γ k∆kα−1 −(q,t)

∆ 
j
j
n
o
2
2
2
2
kβkα-(p,s) k∆kα−1 -(q,t) + 2 kβkα-(p,s) |ei | k∆kα−1 -(q,t) − γ k∆kα−1 -(q,t)
=ei
sup
k∆kα−1 -(q,t)

(
(14)

=

e2i γ−kβkγ2

α-(p,s)

+∞

if γ > kβk2α-(p,s) ,
if γ ≤ kβk2α-(p,s) .

,

where the development uses the duality results developed in Proposition 1. The last equality
is optimize over ∆ for two different cases of λ.
Since optimization over γ is a minimization, the outer player will always select γ that avoids
an infinite value of the game. Then we can write the worst-case expected loss function as,
h
2 i
T
sup
EP Y − X β
(15)
P :Dc (P,Pn )≤δ

)
EPn l (X, Y ; β)
= min
γδ − γ
γ>kβk2α-(p,s)
γ − kβk2α-(p,s)
2
p
√
EPn l (X, Y ; β) + δ kβkα-(p,s) .
=
(

The first equality in (15) is a plug-in from the result in (14). For the second equality, we can
observe the target function is convex and differentiable and as γ → ∞ and γ → kβk2α-(p,s) ,
the value function will be infinity. We can solve this convex optimization problem which
leads to the result above. We further take square root and take minimization over β on both
sides, we proved the claim of the theorem.

A.3. Proof for Optimal Selection of Regularization for Linear Regression.
Proof for Theorem 3. For linear regression with the square loss function, if we apply the
strong duality result for semi-infinity linear programming problem as in Section B of Blanchet
et al. (2016), we can write the scaled RWP function for linear regression as

(16)
nRn (β ∗ ) = sup −ζ T Zn − EPn φ (Xi , Yi , β ∗ , ζ) ,
ζ
√1
n

Pn

and
n
 T  
 T o
2
∗
T
∗T
−1/2
∗T
φ (Xi , Yi , β , ζ) = sup ei ζ ∆ − β ∆ ζ Xi − k∆kα−1 -(q,t) + n
β ∆ ζ ∆
.

where Zn =

i=1 ei Xi

∆

Follow the similar discussion in the proof of Theorem 4 in Blanchet et al. (2016). Applying
Lemma 2 in Blanchet et al. (2016), we can argue that the optimizer ζ can be restrict on a
compact set asymptotically with high probability. We can apply the uniform law of large

16

number estimate as in Lemma 3 of Blanchet et al. (2016) to the second term in (16) and we
obtain
nRn (β ∗ ) = sup{−ζ T Zn − Eφ (X, Y, β, ζ)]} + oP (1).

(17)

ζ

For any fixed X, Y , as n → ∞, we can simplify the contribution of φ (·) inside sup in (17).
This is done by applying the duality result (Hölder-type inequality) in Proposition 1 and
noting that φ (·) becomes quadratic in k∆kα−1 -(q,t) . This results in the simplified expression
io
n
h
2
nRn (β ∗ ) = sup −ζ T Zn − E eζ − (ζ T X)β ∗ α-(p,s) + oP (1).
ζ

Since we can observe that, Zn ⇒ σZ, then as n → ∞ we proved the first argument. For
this step we need to show that the feasible region can be compactified with high probability.
This compactification argument is done using a technique similar to Lemma 2 in Blanchet
et al. (2016).
By the definition of L1 , we can apply Hölder inequality to the first term, and split the
optimization into optimizing over direction kζ 0 kα-(p,s) = 1 and magnitude a ≥ 0. Thus, we
have
io
n
h
2
L1 ≤ 0 0max
max 2aσ kZkα−1 -(q,t) − a2 E eζ 0 − (ζ 0T X)β ∗ α-(p,s) .
ζ :kζ kα-(p,s) =1 a≥0

It is easy to solve the quadratic programming problem in a and we conclude that
σ 2 kZk2α−1 -(q,t)
h
i.
L1 ≤
2
0
0T
∗
minζ 0 :kζ 0 kα-(p,s) =1 E keζ − (ζ X)β kα-(p,s)
For the denominator, we have estimates as follows:
h
i
0
0T
∗ 2
eζ
−
(ζ
X)β
min
E
α-(p,s)
0
0
ζ :kζ kα-(p,s) =1

≥

h

min

ζ 0 :kζ 0 kα-(p,s) =1

≥ V ar(|e|) +

T

∗

E |e| − ζ X kβ kα-(p,s)
min
0

ζ 0 :kζ kα-(p,s) =1



i2

2
kβ ∗ kα-(p,s) E ζ 0T X − E |e| ≥ V ar(|e|).

The first estimate is due to the triangle inequality in Proposition 1, the second estimate follows using Jensen’s inequality, the last inequality is immediate. Combining these inequalities
we conclude
L1 ≤ σ 2 kZk2α−1 -(q,t) /V ar(|e|).

Appendix B. Additional Materials
In this Section, we will provide the proofs for DRO representation and asymptotic result
for logistic regression, which were discussed in Theorem 2 and Theorem 4, in Section B.1
and Section B.2. In addition, we will provide the results under the high dimension setting
for linear regression, where the number of predictors growth with the sample size, as a
generalization of Theorem 3, which we proved in Section B.3.

17

B.1. Proof of DRO for Logistic Regression.
Proof for Theorem 2. By applying strong duality results for semi-infinity linear programming
problem in Blanchet et al. (2016), we can write the worst case expected loss function as,


sup
EP log 1 + exp −Y β T X
P :Dc (P,Pn )≤δ

)
n
o
n

1X
.
sup log 1 + exp −Yi β T u − γ kXi − ukα−1 -(q,t)
= min γδ −
γ≥0
n i=1 u
(

For each i, we can apply Lemma 1 in Shafieezadeh-Abadeh et al. (2015) and the dual norm
result in Proposition 1 to deal with the inner optimization problem. It gives us,
n
o

sup log 1 + exp −Yi β T u − γ kXi − ukα−1 -(q,t)
u


log 1 + exp −Yi β T Xi
if kβkα-(p,s) ≤ γ,
=
∞
if kβkα-(p,s) > γ.
Moreover, since the outer player wishes to minimize, γ will be chosen to satisfy γ ≥ kβkα-(p,s) .
We then conclude
(
)
n
n
o

1X
T
min γδ −
sup log 1 + exp −Yi β u − γ kXi − ukα−1 -(q,t)
γ≥0
n i=1 u
(
)
n

1X
= min
δγ +
log 1 + exp −Yi β T Xi
γ≥kβkα-(p,s)
n i=1
n


1X
=
log 1 + exp −Yi β T Xi + δ kβkα-(p,s) ,
n i=1
where the last equality is obtained by noting that the objective function is continuous and
monotone increasing in γ, thus γ = kβkα-(p,s) is optimal. Hence, we conclude the DRO
formulation for GR-Lasso logistic regression.

B.2. Proof of Optimal Selection of Regularization for Logistic Regression.
Proof of Theorem 4. We can apply strong duality result for semi-infinite linear programming
problem in Section B of Blanchet et al. (2016), and write the scaled RWP function evaluated
at β ∗ in the dual form as,

√
nRn (β ∗ ) = max ζ T Zn − EPn φ (X, Y, β ∗ , ζ) ,
ζ

where Zn =

1
n

Pn

Yi Xi
i 1+exp(Yi X T β ∗ )
i



T
φ (X, Y, β , ζ) = max Y ζ
∗

u

and

X
u
−
T
∗
1 + exp (Y X β ) 1 + exp (Y uT β ∗ )




− kX − ukα−1 -(q,t) .

We proceed as in our proof of Theorem 3 in this paper and also adapting the case ρ = 1 for
Theorem 1 in Blanchet et al. (2016). We can apply Lemma 2 in Blanchet et al. (2016) and
conclude that the optimizer ζ can be taken to lie within a compact set with high probability

18

as n → ∞. We can combine the uniform law of large number estimate as in Lemma 3 of
Blanchet et al. (2016) and obtain

√
nRn (β) = max ζ T Zn − EP φ (X, Y, β ∗ , ζ) + oP (1).
ζ

For the optimization problem defining φ (·), we can apply results in Lemma 5 in Section A.3
of Blanchet et al. (2016), we know, for any choice of ζ̃, if,

T ∗
Id×d − XX T
T y 1 + exp Y X β
ess sup ζ̃
> 1,
(1 + exp (Y X T β ∗ ))2
X,Y
α-(p,s)
h 
i
we have E φ X, Y, β ∗ , ζ̃
= ∞. Since the outer optimization problem is maximization
over ζ, the player will restrict ζ within the set A, where





T
T ∗
Id×d − XX
y 1 + exp Y X β
≤
1
.
A = ζ ∈ Rd : ess sup ζ T


(1 + exp (Y X T β ∗ ))2
X,Y
α-(p,s)
Moreover, it is easy to calculate, if ζ ∈ A, we have E[φ (X, Y, β ∗ , ζ)] = 0, thus we have the
scaled RWP function has the following estimate, as n → ∞
√
nRn (β) = max ζ T Zn + oP (1).
ζ∈A

Letting n → ∞, we obtain the exact asymptotic result.
For the stochastic upper bound, let us recall for the definition of the set A and consider the
following estimate

T ∗
y
1
+
exp
Y
X
β
Id×d − XX T
T
ζ
(1 + exp (Y X T β ∗ ))2
α-(p,s)
≥

Yζ
1 + exp (Y β ∗ T X)

−
α-(p,s)

ζ T Xβ ∗
(1 + exp (Y β ∗ T X))2
∗

≥

α-(p,s)

kXkα−1 -(q,t) kβ kα-(p,s)
1
−
1 + exp (Y β ∗ T X) (1 + exp (Y β ∗ T X)) (1 + exp (−Y β ∗ T X))

!
kζkα-(p,s) .

The first inequality is due to application of triangle inequality in Proposition 1, while the
second estimate follows from Hölder’s inequality and Y ∈ {−1, +1}. Since we assume positive
probability density for the predictor X, we can argue that, if kζkα-(p,s) = (1 − )−2 > 1 and
 > 0 is chosen arbitrarily small, we can conclude from the above estimate that, we have

T ∗
Id×d − XX T
T y 1 + exp Y X β
> 1.
ζ
(1 + exp (Y X T β ∗ ))2
α-(p,s)

19

n
o
Thus, we proved the claim that A ⊂ ζ, kζkα-(p,s) ≤ 1 . The stochastic upper bound is
n
o
derived by replacing A by ζ, kζkα-(p,s) ≤ 1 , i.e.
L3 = sup ζ T Z ≤

sup
kζkα-(p,s) ≤1

ζ∈A

ζ T Z = kZkα−1 -(q,t) ,

where the final estimation is due to dual norm structure in Proposition 1. Since we know,
1
≤ 1, it is easy to argue, V ar(Z̃) − V ar(Z) is positive semidefinite, thus, we know
1+exp(Y X T β)
kZkα−1 -(q,t) is stochastic dominated by L4 := Z̃

α−1 -(q,t)

. Hence, we obtain L3 ≤ L4 .



B.3. Technical Results for Optimal Regularization in GSRL for High Dimensional
Linear Regression. We conclude the section by exploring the behavior of the optimal
distributional uncertainty (in the sense of optimality presented in Section 3) as the dimension
increases. This is an analog of the high-dimension result for SR-Lasso as Theorem 6 in
Blanchet et al. (2016).
Theorem 5 (RWP Function Asymptotic Results for High-dimension). Suppose that as√
√ T
sumptions in Theorem 3 hold and select p = 2, s = 1 let us writeg̃ −1 =
g1 , . . . , gd¯ (so
√
√ T
√
αj = gj ) and g̃ −1/2 = 1/ g1 , . . . , 1/ gd¯ respectively. Moreover, let us define C (n, d)
h
i
d¯ √
√
E
max
g
kX
(G
)k
i
i 2
E kXk d-(2,1)
i=1
¯
√
√
C (n, d) =
=
.
n
n

Assume that largest eigenvalue of Σ is of order o nC (n, d)2 , that β ∗ satisfies a weak sparsity
condition, namely, kβ ∗ k√g̃-(2,1) = o (1/C (n, d)). Then,
nRn (β ∗ ) .D
as n, d → ∞, where Zn := n−1/2

kZn kg̃−1/2 -(2,∞)
V ar (|e|)

,

Pn

i=1 ei Xi .

Proof. For linear regression model with square loss function, the RWP function is defined
as in equation (10). By considering the cost function as in Theorem 1 and applying the
strong duality results in the Appendix of Blanchet et al. (2016), we can write the scaled
RWP function in the dual form as,

nRn (β ∗ ) = sup − ζ T Zn
ζ
n

 √


1 X
−√
sup{ei ζ T ∆ − β ∗ T ∆ ζ T Xi −
n k∆k2g̃−1/2 -(2,∞) + β ∗ T ∆ ζ T ∆ } .
n i=1 ∆

For each i−th inneroptimization
problem, we can apply Hölder inequality in Proposition 1

for the term β ∗ T ∆ ζ T ∆ , we have an upper bound for the scaled RWP function, i.e.

nRn (β ∗ ) ≤ sup − ζ T Zn
ζ
!
n
∗
kβ
k


√
−1 -(p,s) kζkg̃ −1 -(p,s)
1 X
T
g̃
√
sup{ ei ζ − ζ T Xi β ∗ ∆ − n 1 −
k∆k2g̃−1/2 -(q,t) } .
−√
n i=1 ∆
n

20

Since the coefficients for each inner optimization problem is negative and we can get an
upper bound for RWP function if we do not fully optimize the inner optimization problem.
For each i, let us take ∆ to the direction satisfying the Hölder inequality in Property 1 for
 T
the term ei ζ − ζ T Xi β ∗ ∆ and only optimize the magnitude of ∆, for simplicity let us
denote γ = k∆kg̃−1/2 -(q,t) .
We have,

nRn (β ∗ ) ≤ sup − ζ T Zn
ζ
!
n
∗
√
kβ
k

√
−1 -(p,s) kζk g̃-(p,s)
1 X
g̃
√
sup{γ ei ζ − ζ T Xi β ∗ √g-(p,s) − n 1 −
γ 2} .
−√
n i=1 γ
n
For each inner optimization problem it is of quadratic form in γ, especially, when n is large
the coefficients for the second order term will be negative, thus, as n → ∞, we can solve the
inner optimization problem and obtain,
nRn (β ∗ )
n


1X
1

≤ sup − ζ Zn − 
e i ζ − ζ T Xi β∗
ζ
4 1 − kβ ∗ kg−1 -(p,s) kζkg̃−1 -(p,s) n−1/2 n i=1


T

2
g̃ −1 -(p,s)

n

= sup



sup

a≥0 ζ:kζk√g̃-(p,s) =1


a2
1X

− aζ Zn − 
ei ζ − ζ T Xi β ∗
4 1 − kβ ∗ kg̃−1 -(p,s) an−1/2 n i=1
T

2
g̃ −1 -(p,s)

The equality above is due to changing to polar coordinate for the ball under g̃ −1 -(p, s) norm.
For the first term, ζ T Zn , when kζkg̃−1 -(p,s) = 1, we can apply Hölder inequality again, i.e.
ζ T Zn ≤ kZn kg̃−1/2 -(q,t) . Then, only the second term in the previous display involves the
direction of ζ, thus we can have

nRn (β ∗ ) ≤ sup a kZn kg̃−1/2 -(q,t)
a≥0

a2

− 
4 1 − kβ ∗ k√g-(p,s)

n


1X

inf
e i ζ − ζ T Xi β ∗
an−1/2 ζ:kζkg̃−1 -(p,s) =1 n i=1

2
√

g-(p,s)

.

By the weak sparsity assumption, we have kβ ∗ kg̃−1 -(p,s) n−1/2 → 0 as n → ∞, the supremum
over a is attained at
2 kZn kg̃−1/2 -(q,t)
+ o(1),
a∗ =
P
inf ζ:kζk√g̃-(p,s) =1 n1 ni=1 kei ζ − (ζ T Xi ) β ∗ k2g−1 -(p,s)
as n → ∞. Therefore, we have the upper bound estimator for the scaled RWP function,
(18)

∗

nRn (β ) ≤

inf ζ:kζkg̃−1 -(p,s) =1

kZn k2g̃−1/2 -(q,t)
+ op (1).
Pn
2
1
T
i=1 kei ζ − (ζ Xi ) β∗ kg −1 -(p,s)
n

To get the final result, we try to find a lower bound for the infimum in the denominator. For
the objective function in the denominator, since we optimize on the surface kζkg̃−1 -(p,s) = 1,

.

21

and due to the triangle inequality analysis in Proposition 1, we have
n

1X
2
ei ζ − ζ T Xi β ∗ g̃−1 -(p,s)
n i=1
n
2
1 X
T
∗ √
≥
|ei | kζkg̃−1 -(p,s) − ζ Xi kβ k g̃-(p,s)
n i=1
n
n
n
1X T
1X T
1X 2
2
|ei | + kβ ∗ k2g̃−1 -(p,s)
ζ Xi − 2 kβ ∗ kg̃−1 -(p,s) E [|ei |]
ζ Xi − n (ζ) ,
n i=1
n i=1
n i=1
P
where n (ζ) = 2 kβ ∗ kg̃−1 -(p,s) n1 ni=1 (|ei | − E [|ei |]). Let us denote the pseudo error to be
ẽi = |ei | − E [|ei |], which has mean zero and V ar [ẽi ] ≤ V ar [ei ]. Since ei is independent of
Xi we have that


E ẽi ζ T Xi = 0,


V ar ẽi ζ T Xi = V ar [ẽi ] ζ T Σζ ≤ V ar [ei ] ζ T Σζ.

=

By our assumptions on the eigenstructure of Σ, i.e. λmax (Σ) = o (nC(n, d)2 ), for the case
p = 2 and s = 1, we have

sup
ζ T Σζ ≤
sup
λmax (Σ) kζk2 ≤ λmax (Σ) = o nC(n, d)2 .
ζ:kζkg̃−1 -(2,1) =1

ζ:kζkg̃−1 -(2,1) =1

P
Then, we have the variance of n1 ni=1 ζ T Xi is of order o (C(n, d)2 ) uniformly on kζkg̃−1 -(p,s) =
1. Combining this estimate with the weak sparsity assumption that we have imposed, we
have
n (ζ) = op (1).
√
Since the estimate is uniform over kζk g̃-(2,1) = 1, we have that for n sufficiently large,
n


1X
ei ζ − ζ T Xi β ∗
n i=1

2
g̃ −1 -(2,1)

!2
n
n
X
1
1X 2
|ei | − (E [|ei |])2 +
inf
kβ ∗ kg̃−1 -(2,1)
ζ T Xi − E [|ei |] + op (1)
≥
ζ:kζkg̃−1 -(2,1) =1
n i=1
n i=1
≥ V arn [|ei |] + op (1).
Combining the above estimate and equation (18), when p = q = 2, s = 1 and t = ∞, we
have that
kZn k2g̃−1/2 -(2,∞)
∗
nRn (β ) ≤
+ op (1),
V ar [|e|]
as n → ∞.


