Provably tuning the ElasticNet across instances

Abstract

An important unresolved challenge in the theory of regularization is to set the regularization coefficients of popular techniques like the ElasticNet with general provable guarantees. We consider the problem of tuning the regularization parameters of Ridge regression, LASSO, and the ElasticNet across multiple problem instances, a setting that encompasses both cross-validation and multi-task hyperparameter optimization. We obtain a novel structural result for the ElasticNet which characterizes the loss as a function of the tuning parameters as a piecewise-rational function with algebraic boundaries. We use this to bound the structural complexity of the regularized loss functions and show generalization guarantees for tuning the ElasticNet regression coefficients in the statistical setting. We also consider the more challenging online learning setting, where we show vanishing average expected regret relative to the optimal parameter pair. We further extend our results to tuning classification algorithms obtained by thresholding regression fits regularized by Ridge, LASSO, or ElasticNet. Our results are the first general learning-theoretic guarantees for this important class of problems that avoid strong assumptions on the data distribution. Furthermore, our guarantees hold for both validation and popular information criterion objectives.

1 Introduction

Ridge regression [HK70, TA77], LASSO [Tib96], and their generalization the ElasticNet [HTF09] are among the most popular algorithms in machine learning and statistics, with applications to linear classification, regression, data analysis, and feature selection [Cha92, ZY07, HTF09, DW18, FDSC${}^{+}$19]. Given a supervised dataset $(X,y)\in\mathbb{R}^{m\times p}\times\mathbb{R}^{m}$ with $m$ datapoints and $p$ features, these algorithms compute the linear predictor

$$\hat{\beta}_{\lambda_{1},\lambda_{2}}^{(X,y)}=\operatorname*{argmin}_{\beta\in\mathbb{R}^{p}}\|y-X\beta\|_{2}^{2}+\lambda_{1}\|\beta\|_{1}+\lambda_{2}\|\beta\|_{2}^{2}.$$ (1)

Here $\lambda_{1},\lambda_{2}\geq 0$ are regularization coefficients constraining the $\ell_{1}$ and $\ell_{2}$ norms, respectively, of the model $\beta$. For general $\lambda_{1}$ and $\lambda_{2}$ the above algorithm is the ElasticNet, while setting $\lambda_{1}=0$ recovers Ridge and setting $\lambda_{2}=0$ recovers LASSO.

These coefficients play a crucial role across fields: in machine learning controlling the norm of $\beta$ implies provable generalization guarantees and prevent over-fitting in practice [MRT12], in data analysis their combined use yields parsimonious and interpretable models [HTF09], and in Bayesian statistics they correspond to imposing specific priors on $\beta$ [Mur12, LL10]. In practice, $\lambda_{2}$ regularizes $\beta$ by uniformly shrinking all coefficients, while $\lambda_{1}$ encourages the model vector to be sparse. This means that while they do yield learning-theoretic and statistical benefits, setting them to be too high will cause models to under-fit the data. The question of how to set the regularization coefficients becomes even more unclear in the case of the ElasticNet, as one must juggle trade-offs between sparsity, feature correlation, and bias when setting both $\lambda_{1}$ and $\lambda_{2}$ simultaneously. As a result, there has been intense empirical and theoretical effort devoted to automatically tuning these parameters. Yet the state-of-the-art is quite unsatisfactory: proposed work consists of either heuristics without formal guarantees [Gib81, KKM15], approaches that optimize over a finite grid or random set instead of the full continuous domain [CLW16], or analyses that involve very strong theoretical assumptions [Zha09].

In this work, we study a variant on the above well-established and intensely studied formulation. The key distinction is that instead of a single dataset $(X,y)$, we consider a collection of datasets or instances of the same underlying regression problem $(X^{(i)},y^{(i)})$ and would like to learn a pair $(\lambda_{1},\lambda_{2})$ that selects a model in equation (1) that has low loss on a validation dataset. This can be useful to model practical settings, for example where new supervised data is obtained several times or where the set of features may change frequently [DP14]. We do not require all examples across datasets to be i.i.d. draws from the same data distribution, and can capture more general data generation scenarios like cross-validation and multi-task learning [ZY21]. Despite these advantages, we remark that our problem formulation is quite different from the standard single dataset setting, where all examples in the dataset are typically assumed to be drawn independently from the same distribution. Our formulation treats the selection of regularization coefficients as data-driven algorithm design, which is often used to study combinatorial problems [GR17, Bal20].

Our main contribution is a new structural result for the ElasticNet Regression problem, which implies generalization guarantees for selecting ElasticNet Regression coefficients in the multiple-instance setting. In particular, Ridge and LASSO regressions are special cases. We extend our results to obtain low regret in the online learning setting, and to tuning related linear classification algorithms. In summary, we make the following key contributions:

• We formulate the problem of tuning the ElasticNet as a question of learning $\lambda_{1}$ and $\lambda_{2}$ simultaneously across multiple problem instances, either generated statistically or coming online. Our formulation captures relevant settings like cross-validation and multi-task learning.

• We provide a novel structural result (Theorem 2.2) that characterizes the loss of the ElasticNet fit. We show that the hyperparameter space can be partitioned by polynomial curves of bounded degrees into pieces where the loss is a bivariate rational function. The result holds for both the usual ElasticNet validation objective and when it is augmented with information criteria like the AIC or BIC.

• An important consequence of our structural result is a bound on the pseudo-dimension (Definition 5) for the loss function class, which yields strong generalization bounds for tuning $\lambda_{1}$ and $\lambda_{2}$ simultaneously in the statistical learning setting (Theorem 3.2). Informally, for ElasticNet regression problems with at most $p$ parameters, for any problem distribution $\mathcal{D}$, we show that $O\left(\frac{1}{\epsilon^{2}}(p^{2}\log\frac{1}{\epsilon}+\log\frac{1}{\delta})\right)$ problems (datasets) are sufficient to learn an $\epsilon$-approximation to the best $(\lambda_{1},\lambda_{2})$, with probability at least $1-\delta$.

• In the online setting, we show under very mild data assumptions—much weaker than prior work—that the problem satisfies a dispersion condition [BDV18, BDS20]. As a result we can tune all parameters across a sequence of instances appearing online and obtain vanishing regret relative to the optimal parameter in hindsight over the sequence (Theorem 3.3) at the rate $\tilde{O}(1/\sqrt{T})$ wrt the length $T$ of the sequence.

• We show how to extend our results to regularized classifiers that perform thresholding on Ridge, LASSO or ElasticNet regression estimates, again providing strong generalization and online learning guarantees (Theorems 4.1, 4.2).

We include a couple of remarks to emphasize the generality and significance of our results. First, in our multiple-instance formulation the different problem instances need not have the same number of examples, or even the same set of features. This allows us to handle practical scenarios where the set of features changes across datasets, and we can learn parameters that work well on average across multiple different but related regression tasks. Second, by generating problem instances iid from a fixed (training + validation) dataset, we can obtain iterations (training/validation splits) of popular cross-validation techniques (including the popular leave-one-out and Monte Carlo CV) and our result implies that $\tilde{O}(p^{2}/\epsilon^{2})$ iterations are enough to determine an ElasticNet parameter $\hat{\lambda}$ with loss within $\epsilon$ (with high probability) of the optimal parameter $\lambda^{*}$ over the distribution induced by the cross-validation splits.

Key challenges and insights. A major challenge in learning the ElasticNet parameters is that the variation of the solution path as a function of $\lambda_{2}$ is hard to characterize. Indeed the original ElasticNet paper [ZH05] suggests using the heuristic of grid search to learn a good $\lambda_{2}$, even though $\lambda_{1}$ may be exactly optimized by computing full solution paths (for each $\lambda_{2}$). We approach this indirectly by utilizing a characterization of the LASSO solution by [Tib13], which is based on the KKT (Karush–Kuhn–Tucker) optimality conditions, to arrive at a precise piecewise structure for the problem. In more detail, we use these conditions to come up with a set of algebraic curves (polynomial equations in $\lambda_{1}$ and $\lambda_{2}$) of bounded degrees, such that the set of possible discontinuities lie along these curves, and the loss function behaves well (a bounded-degree rational function) in each piece of the partition of the parameter domain induced by these curves. This characterization is crucial in establishing a bound on the structural complexity needed to provide strong generalization guarantees. We further show additional structure on these algebraic curves that (roughly speaking) imply that the curves do not concentrate in any region of the domain, allowing us to use the powerful recipe of [BDP20] for online learning.

1.1 Related work

Model selection for Ridge regression, LASSO or ElasticNet is typically done by selecting the regularization parameter $\lambda$ that works well for given data, although some parameter-free techniques for variable selection have been recently proposed [LM15]. Choosing ‘optimal’ parameters for tuning the regularization has been a subject of extensive theoretical and applied research. Much of this effort is heuristic [Gib81, KKM15] or focused on developing tuning objectives beyond validation accuracy like AIC or BIC [Aka74, Sch78] without providing procedures for provably optimizing them. The standard approach given a tuning objective is to optimize it over a grid or random set of parameters, for which there are guarantees [CLW16], but this does not ensure optimality over the entire continuous tuning domain, especially since objectives such as 0-1 validation error or information criteria can have many discontinuities. Selecting a grid that is too fine or too coarse can result in either very inefficient or highly inaccurate estimates (respectively) for good parameters. Other guarantees make strong assumptions on the data distribution such as sub-Gaussian noise [Zha09, CLC21] or depend on unknown parameters that are hard to quantify in practice [FL10]. Recent work has shown asymptotic consistency of cross-validation for ridge regression, even in the limiting case $\lambda_{2}\rightarrow 0$ which is particularly interesting for the overparameterized regime [HMRT22, PWRT21]. A successful line of work has focused on efficiently obtaining models for different values of $\lambda_{1}$ using regularization paths [EHJT04], but the guarantees are computational rather than learning-theoretic or statistical. In contrast, we provide principled approaches that guarantee near-optimality of selected parameters with high confidence over the entire continuous domain of parameters.

Data-driven algorithm design has proved successful for tuning parameters for a variety of combinatorial problems like clustering, integer programming, auction design and graph-based learning [BDL19, BPSV21, BSV16, BS21]. We provide an application of these techniques to parameter tuning in a problem that is not inherently combinatorial by revealing a novel discrete structure. We identify the underlying piecewise structure of the ElasticNet loss function which is extremely effective in establishing learning-theoretic guarantees [BDD${}^{+}$21]. To exploit this piecewise structure, we analyze the learning-theoretic complexity of rational algebraic function classes and infer generalization guarantees. Follow-up work [BNS23] improves on our generalzation guarantees and extends the results to regularized logistic regression. We also employ and extend general tools and techniques for online data-driven learning from [BDP20, BS21] to rational functions in order to prove our online learning guarantees for regularization coefficient tuning.

2 Preliminaries and a Key Structural Result

Given data $(X,y)$ with $X\in\mathbb{R}^{m\times p}$ and $y\in\mathbb{R}^{m}$, consisting of $m$ labeled examples with $p$ features, we seek estimators $\beta\in\mathbb{R}^{p}$ which minimize the regularized loss. Popular regularization methods like LASSO and ElasticNet can be expressed as computing the solution of an optimization problem given by

$$\hat{\beta}_{\lambda,f}^{(X,y)}\in\operatorname*{argmin}_{\beta\in\mathbb{R}^{p}}\left\lVert y-X\beta\right\rVert_{2}^{2}+\langle\lambda,f(\beta)\rangle,$$

where $f:\mathbb{R}^{p}\rightarrow\mathbb{R}_{\geq 0}^{d}$ gives the regularization penalty for estimator $\beta$, $\lambda\in\mathbb{R}_{\geq 0}^{d}$ is the regularization parameter, and $d$ is the number of regularization parameters. $d=1$ for Ridge and LASSO, and $d=2$ for the ElasticNet. Setting $f=f_{2}$ with $f_{2}(\beta)=\left\lVert\beta\right\rVert_{2}^{2}$ yields Ridge regression, and setting $f(\beta)=f_{1}(\beta):=\left\lVert\beta\right\rVert_{1}$ corresponds to LASSO. Also using $f_{\text{EN}}(\beta):=(f_{1}(\beta),f_{2}(\beta))$ gives the ElasticNet with regularization parameter $\lambda=(\lambda_{1},\lambda_{2})$. Note that we use the same $\lambda$ (with some notational overloading) to denote the regularization parameters for ridge, LASSO, or ElasticNet. We write $\hat{\beta}_{\lambda,f}^{(X,y)}$ as simply $\hat{\beta}_{\lambda,f}$ when the dataset $(X,y)$ is clear from context. On any instance $x\in\mathbb{R}^{p}$ from the feature space, the prediction of the regularized estimator is given by the dot product $\langle x,\hat{\beta}_{\lambda,f}\rangle$. The average squared loss over a dataset $(X^{\prime},y^{\prime})$ with $X^{\prime}\in\mathbb{R}^{m^{\prime}\times p}$ and $y^{\prime}\in\mathbb{R}^{m^{\prime}}$ is given by

$$l_{r}(\hat{\beta}_{\lambda,f},(X^{\prime},y^{\prime}))=\frac{1}{m^{\prime}}\left\lVert y^{\prime}-X^{\prime}\hat{\beta}_{\lambda,f}\right\rVert_{2}^{2}.$$

By setting $(X^{\prime},y^{\prime})$ to be the training data $(X,y)$, we get the training loss $l_{r}(\hat{\beta}_{\lambda,f},(X,y))$. We use $(X_{\text{val}},y_{\text{val}})$ to denote a validation split.

Distributional and Online Settings. In the distributional or statistical setting, we receive a collection of $n$ instances of the regression problem

$$P^{(i)}=(X^{(i)},y^{(i)},X_{\text{val}}^{(i)},y_{\text{val}}^{(i)})\in\mathcal{R}_{m_{i},p_{i},m_{i}^{\prime}}:=\mathbb{R}^{m_{i}\times p_{i}}\times\mathbb{R}^{m_{i}}\times\mathbb{R}^{m^{\prime}_{i}\times p_{i}}\times\mathbb{R}^{m^{\prime}_{i}},$$

for $i\in[n]$ generated i.i.d. from some problem distribution $\mathcal{D}$. The problems are in the problem space given by $\Pi_{m,p}=\bigcup_{m_{1}\geq 0,m_{2}\leq m,p_{1}\leq p}\mathcal{R}_{m_{1},p_{1},m_{2}}$ (note that the problem distribution $\mathcal{D}$ is over $\Pi_{m,p}$). On any given instance $P^{(i)}$ the loss is given by the squared loss on the validation set, $\ell_{\text{EN}}(\lambda,P^{(i)})=l_{r}(\hat{\beta}_{\lambda,f_{\text{EN}}}^{(X^{(i)},y^{(i)})},(X_{\text{val}}^{(i)},y_{\text{val}}^{(i)}))$. On the other hand, in the online setting, we receive a sequence of $T$ instances of the ElasticNet regression problem $P^{(i)}=(X^{(i)},y^{(i)},X_{\text{val}}^{(i)},y_{\text{val}}^{(i)})\in\Pi_{m,p}$ for $i\in[T]$ online. On any given instance $P^{(i)}$, the online learner is required to select the regularization parameter $\lambda^{(i)}$ without observing $y_{\text{val}}^{(i)}$, and experiences loss given by $\ell(\lambda^{(i)},P^{(i)})=l_{c}(\hat{\beta}_{\lambda^{(i)},f_{EN}}^{(X^{(i)},y^{(i)})},(X_{\text{val}}^{(i)},y_{\text{val}}^{(i)}))$. The goal is to minimize the regret w.r.t. choosing the best fixed parameter in hindsight for the same problem sequence, i.e.

$$R_{T}=\sum_{i=1}^{T}\ell(\lambda^{(i)},P^{(i)})-\min_{\lambda}\sum_{i=1}^{T}\ell(\lambda,P^{(i)}).$$

We also define average regret as $\frac{1}{T}R_{T}$ and expected regret as $\mathbb{E}[R_{T}]$ where the expectation is over both the randomness of the loss functions and any random coins used by the online algorithm.

Given a class of regularization algorithms $\mathcal{A}$ parameterized by regularization parameter $\lambda$ over a set of problem instances $\mathcal{X}$, and given loss function $\ell:\mathcal{A}\times\mathcal{X}\rightarrow\mathbb{R}$ which measures the loss of any algorithm in $\mathcal{A}$ on any fixed problem instance, consider the set of functions $\mathcal{H}_{\mathcal{A}}=\{\ell(A,\cdot)\mid A\in\mathcal{A}\}$. For example, for the ElasticNet we have $\ell_{\text{EN}}(\lambda,P)=l_{r}(\hat{\beta}^{(X_{P},y_{P})}_{\lambda,f_{\text{EN}}},(X^{\prime}_{P},y^{\prime}_{P}))$, where $(X_{P},y_{P})$ and $(X^{\prime}_{P},y^{\prime}_{P})$ are the training and validation sets associated with problem $P\in\mathcal{X}$ respectively. Bounding the pseudo-dimension of $\mathcal{H}_{\mathcal{A}}$ gives a bound on the sample complexity for uniform convergence guarantees, i.e. a bound on the sample size $n$ for which the algorithm $\hat{A}_{S}\in\mathcal{A}$ which minimizes the average loss on any sample $S$ of size $n$ drawn i.i.d. from any problem distribution $\mathcal{D}$ is guaranteed to be near-optimal with high probability [Dud67]. See Appendix A for the relevant classic definitions and results. Define the dual class $\mathcal{H}^{*}$ of a set of real-valued functions $\mathcal{H}\subseteq 2^{\mathcal{X}}$ as $\mathcal{H}^{*}=\{h^{*}_{x}:\mathcal{H}\rightarrow\mathbb{R}\mid x\in\mathcal{X}\}$ where $h^{*}_{x}(h)=h(x)$. In the context of regression problems $\mathcal{X}$, for each fixed problem instance $x\in\mathcal{X}$ there is a dual function $h^{*}_{x}$ that computes the loss $\ell(A,x)$ for any (primal) function $h_{A}=\ell(A,\cdot)\in\mathcal{H}_{\mathcal{A}}$. For a function class $\mathcal{H}$, showing that dual class $\mathcal{H}^{*}$ is piecewise-structured in the sense of Definition 1 and bounding the complexity of the duals of boundary and piece functions of $\mathcal{H}^{*}$ are useful to understand the learnability of $\mathcal{H}$ [BDD${}^{+}$21].

Definition 1 (Piecewise structured functions, [BDD${}^{+}$21]).

A function class $H\subseteq\mathbb{R}^{\mathcal{X}}$ that maps a domain $\mathcal{X}$ to $\mathbb{R}$ is $(F,G,k)$-piecewise decomposable for a class $G\subseteq\{0,1\}^{\mathcal{X}}$ of boundary functions and a class $F\subseteq\mathbb{R}^{\mathcal{X}}$ of piece functions if the following holds: for every $h\in H$, there are $k$ boundary functions $g_{1},\dots,g_{k}\in G$ and a piece function $f_{\mathbf{b}}\in F$ for each bit vector $\mathbf{b}\in\{0,1\}^{k}$ such that for all $x\in\mathcal{X}$, $h(x)=f_{\mathbf{b}_{x}}(x)$ where $\mathbf{b}_{x}=(g_{1}(x),\dots,g_{k}(x))\in\{0,1\}^{k}$.

Intuitively, a real-valued function is piecewise-structured if the domain can be divided into pieces by a finite number of boundary functions (say linear or polynomial thresholds) and the function value over each piece is easy to characterize (e.g. constant, linear, polynomial). To state and understand our structural insights into the ElasticNet problem we will also need the definition of equicorrelation sets, the subset of features with maximum absolute correlation for any fixed $\lambda_{1}$, useful for characterizing LASSO/ElasticNet solutions. For any subset $\mathcal{E}\subseteq[p]$ of the features, we define $X_{\mathcal{E}}=\begin{pmatrix}\dots X_{*i}\dots\end{pmatrix}_{i\in\mathcal{E}}$ as the $m\times|\mathcal{E}|$ matrix of columns $X_{*i}$ of $X$ corresponding to indices $i\in\mathcal{E}$. Similarly $\beta_{\mathcal{E}}\in\mathbb{R}^{|\mathcal{E}|}$ is the subset of estimators in $\beta$ corresponding to indices in $\mathcal{E}$. We will assume all the feature matrixes $X$ (for training datasets) are in general position (Definition 6).

Definition 2 (Equicorrelation sets, [Tib13]).

Let $\beta^{*}\in\operatorname*{argmin}_{\beta\in\mathbb{R}^{p}}\left\lVert y-X\beta\right\rVert_{2}^{2}+\lambda_{1}||\beta||_{1}$. The equicorrelation set corresponding to $\beta^{*}$, $\mathcal{E}=\{j\in[p]\mid|\bm{x}_{j}^{T}(y-X\beta^{*})|=\lambda_{1}\}$, is simply the set of covariates with maximum absolute correlation. We also define the equicorrelation sign vector for $\beta^{*}$ as $s=\textup{{sign}}(X_{\mathcal{E}}^{T}(y-X\beta^{*}))\in\{\pm 1\}$.

Here $\textup{{sign}}(x)=1$ if $x\geq 0$, and $\textup{{sign}}(x)=-1$ otherwise. Consider the class of algorithms consisting of ElasticNet regressors for different values of $\lambda=(\lambda_{1},\lambda_{2})\in(0,\infty)\times(0,\infty)$. We assume $\lambda_{1}>0$ for technical simplicity (cf. [Tib13]). We seek to solve problems of the form $P=(X,y,X_{\text{val}},y_{\text{val}})\in\Pi_{m,p}$, where $(X,y)$ is the training set, $(X_{\text{val}},y_{\text{val}})$ is the validation set with the same set of features, and $m,p$ are upper bounds on the number of examples and features respectively in any dataset. Let $\mathcal{H}_{\text{EN}}=\{\ell_{\text{EN}}(\lambda,\cdot)\mid\lambda\in(0,\infty)\times(0,\infty)\}$ denote the set of loss functions for the class of algorithms consisting of ElasticNet regressors for different values of $\lambda\in\mathbb{R}^{+}\times\mathbb{R}^{+}$. Additionally, we will consider information criterion based loss functions, $\ell_{\text{EN}}^{\textsf{AIC}}(\lambda,P)=\ell_{\text{EN}}(\lambda,P)+2||\hat{\beta}_{\lambda,f_{\text{EN}}}^{(X,y)}||_{0}$ and $\ell_{\text{EN}}^{\textsf{BIC}}(\lambda,P)=\ell_{\text{EN}}(\lambda,P)+2||\hat{\beta}_{\lambda,f_{\text{EN}}}^{(X,y)}||_{0}\log m$ [Aka74, Sch78]. Let $\mathcal{H}_{\text{EN}}^{\textsf{AIC}}$ and $\mathcal{H}_{\text{EN}}^{\textsf{BIC}}$ denote the corresponding sets of loss functions. These criteria are popularly used to compute the squared loss on the training set, to give alternatives to cross-validation. We do not make any assumption on the relation between training and validation sets in our formulation, so our analysis can capture these settings as well.

2.1 Piecewise structure of the ElasticNet loss

We will now establish a piecewise structure of the dual class loss functions (Definition 1). A key observation is that if the signed equicorrelation set $(\mathcal{E},s)$ (i.e. a subset of features $\mathcal{E}\subseteq[p]$ with the same maximum absolute correlation, assigned a fixed sign pattern $\{-1,+1\}^{|\mathcal{E}|}$, see Definition 2) is fixed, then the ElasticNet coefficients may be characterized (Lemma C.1) and the loss is a fixed rational polynomial piece function of the parameters $\lambda_{1},\lambda_{2}$. We then show the existence of a set of boundary function curves $\mathcal{G}$, such that any region of the parameter space located on a fixed side of all the curves (more formally, for a fixed sign pattern in Definition 1) in $\mathcal{G}$ has the same signed equicorrelation set. The boundary functions are a collection of possible curves at which a covariate may enter or leave the set $\mathcal{E}$ and correspond to polynomial thresholds. We make repeated use of the following lemma which provides useful properties of the piece functions as well the the boundary functions of the dual class loss functions.

Figure 1: An illustration of the piecewise structure of the ElasticNet loss, as a function of the regularization parameters, for a fixed problem instance. Pieces are regions where some bounded degree polynomials ($r_{1},r_{2}$) have a fixed sign pattern (one of $\pm 1,\pm 1$), and in each piece the loss is a fixed (rational) function.

Lemma 2.1.

Let $A$ be an $r\times s$ matrix. Consider the matrix $B(\lambda)=(A^{T}A+\lambda I_{s})^{-1}$ and $\lambda>0$.

• Each entry of $B(\lambda)$ is a rational polynomial $P_{ij}(\lambda)/Q(\lambda)$ for $i,j\in[s]$ with each $P_{ij}$ of degree at most $s-1$, and $Q$ of degree $s$.

• Further, for $i=j$, $P_{ij}$ has degree $s-1$ and leading coefficient 1, and for $i\neq j$ $P_{ij}$ has degree at most $s-2$. Also, $Q(\lambda)$ has leading coefficient $1$.

The proof is straightforward and deferred to Appendix C. We will now formally state and prove our key structural result which is needed to establish our generalization and online regret guarantees in Section 3.

Theorem 2.2.

Let $\mathcal{L}$ be a set of functions $\{l_{\lambda}:\Pi_{m,p}\rightarrow\mathbb{R}_{\geq 0}\mid\lambda\in\mathbb{R}^{+}\times\mathbb{R}^{+}$ that map a regression problem instance $P\in\Pi_{m,p}$ to the validation loss $\ell_{\text{EN}}(\lambda,P)$ of ElasticNet trained with regularization parameter $\lambda=(\lambda_{1},\lambda_{2})$. The dual class $\mathcal{L}^{*}$ is $(\mathcal{F},\mathcal{G},p3^{p})$-piecewise decomposable, with $\mathcal{F}=\{f_{q}:\mathcal{L}\rightarrow\mathbb{R}\}$ consisting of rational polynomial functions $f_{q_{1},q_{2}}:l_{\lambda}\mapsto\frac{q_{1}(\lambda_{1},\lambda_{2})}{q_{2}(\lambda_{2})}$, where $q_{1},q_{2}$ have degrees at most $2p$, and $\mathcal{G}=\{g_{r}:\mathcal{L}\rightarrow\{0,1\}\}$ consists of polynomial threshold functions $g_{r}:u_{\lambda}\mapsto\mathbb{I}\{r(\lambda_{1},\lambda_{2})<0\}$, where $r$ is a polynomial of degree 1 in $\lambda_{1}$ and at most $p$ in $\lambda_{2}$.

Proof.

Let $P=(X,y,X_{\text{val}},y_{\text{val}})\in\Pi_{m,p}$ be a regression problem instance. By using the standard reduction to LASSO [ZH05] and well-known characterization of the LASSO solution in terms of equicorrelation sets, we can characterize the solution $\hat{\beta}_{\lambda,f_{EN}}$ of the Elastic Net as follows (Lemma C.1):

$$\hat{\beta}_{\lambda,f_{EN}}=(X_{\mathcal{E}}^{T}X_{\mathcal{E}}+\lambda_{2}I_{|\mathcal{E}|})^{-1}X_{\mathcal{E}}^{T}y-\lambda_{1}(X_{\mathcal{E}}^{T}X_{\mathcal{E}}+\lambda_{2}I_{|\mathcal{E}|})^{-1}s,$$

for some $\mathcal{E}\in[p]$ and $s\in\{-1,1\}^{p}$. Thus for any $\lambda=(\lambda_{1},\lambda_{2})$, the prediction $\hat{y}$ on any validation example with features $\bm{x}\in\mathbb{R}^{p}$ satisfies (for some $\mathcal{E},s\in 2^{[p]}\times\{-1,1\}^{p}$)

$$\hat{y}_{j}=\bm{x}\hat{\beta}_{\lambda,f_{EN}}=\bm{x}(X_{\mathcal{E}}^{T}X_{\mathcal{E}}+\lambda_{2}I_{|\mathcal{E}|})^{-1}X_{\mathcal{E}}^{T}y-\lambda_{1}\bm{x}(X_{\mathcal{E}}^{T}X_{\mathcal{E}}+\lambda_{2}I_{|\mathcal{E}|})^{-1}s.$$

For any subset $R\subseteq\mathbb{R}^{2}$, if the signed equicorrelation set $(\mathcal{E},s)$ is fixed over $R$, then the above observation, together with Lemma 2.1 implies that the loss function $\ell_{\text{EN}}(\lambda,P)$ is a rational function of the form $\frac{q_{1}(\lambda_{1},\lambda_{2})}{q_{2}(\lambda_{2})}$, where $q_{1}$ is a bivariate polynomial with degree at most $2|\mathcal{E}|$ and $q_{2}$ is univariate with degree $2|\mathcal{E}|$.

To show the piecewise structure, we need to demonstrate a set boundary functions $\mathcal{G}=\{g_{1},\dots,g_{k}\}$ such that for any sign pattern $\mathbf{b}\in\{0,1\}^{k}$, the signed equicorrelation set $(\mathcal{E},s)$ for the region with sign pattern $\mathbf{b}$ is fixed. To this end, based on the observation above, we will consider the conditions (on $\lambda$) under which a covariate may enter or leave the equicorrelation set. We will show that this can happen only at one of a finite number of algebraic curves (with bounded degrees).

Condition for joining $\mathcal{E}$. Fix $\mathcal{E},s$. Also fix $j\notin\mathcal{E}$. If covariate $j$ enters the equicorrelation set, the KKT conditions (Lemma B.1) applied to the LASSO problem corresponding to the ElasticNet (Lemma C.1) imply

$$(\bm{x}_{j}^{*})^{T}(y^{*}-X_{\mathcal{E}}^{*}(c_{1}-c_{2}\lambda_{1}^{*}))=\pm\lambda_{1}^{*},$$

where $c_{1}=({X^{*}_{\mathcal{E}}}^{T}X^{*}_{\mathcal{E}})^{-1}{X^{*}_{\mathcal{E}}}^{T}y^{*}$, $c_{2}=({X^{*}_{\mathcal{E}}}^{T}X^{*}_{\mathcal{E}})^{-1}s$, $X^{*}=\frac{1}{\sqrt{1+\lambda_{2}}}\begin{pmatrix}X\\ \sqrt{\lambda_{2}}I_{p}\end{pmatrix}$, $y^{*}=\begin{pmatrix}y\\ 0\end{pmatrix}$, and $\lambda_{1}^{*}=\frac{\lambda_{1}}{\sqrt{1+\lambda_{2}}}$. Rearranging, and simplifying, we get

$$\lambda_{1}^{*} =\frac{(\bm{x}_{j}^{*})^{T}X^{*}_{\mathcal{E}}({X^{*}_{\mathcal{E}}}^{T}X^{*}_{\mathcal{E}})^{-1}(X^{*}_{\mathcal{E}})^{T}y^{*}-(\bm{x}^{*}_{j})^{T}y^{*}}{(\bm{x}_{j}^{*})^{T}X^{*}_{\mathcal{E}}({X^{*}_{\mathcal{E}}}^{T}X^{*}_{\mathcal{E}})^{-1}s\pm 1},\text{or}$$

$$\lambda_{1} =\frac{\bm{x}_{j}^{T}X_{\mathcal{E}}({X_{\mathcal{E}}}^{T}X_{\mathcal{E}}+\lambda_{2}I_{|\mathcal{E}|})^{-1}{X_{\mathcal{E}}}^{T}y-\bm{x}_{j}^{T}y}{\bm{x}_{j}^{T}X_{\mathcal{E}}({X_{\mathcal{E}}}^{T}X_{\mathcal{E}}+\lambda_{2}I_{|\mathcal{E}|})^{-1}s\pm 1}.$$

Note that the terms $(\bm{x}_{j}^{*})^{T}X^{*}_{\mathcal{E}}=\bm{x}_{j}^{T}X_{\mathcal{E}}$, $(X^{*}_{\mathcal{E}})^{T}y^{*}=X_{\mathcal{E}}^{T}y$, and $(\bm{x}^{*}_{j})^{T}y^{*}=\bm{x}_{j}^{T}y$ do not depend on $\lambda_{1}$ or $\lambda_{2}$ (the $\lambda_{2}$ terms are zeroed out since $j\notin\mathcal{E}$). Moreover, $({X^{*}_{\mathcal{E}}}^{T}X^{*}_{\mathcal{E}})^{-1}=({X_{\mathcal{E}}}^{T}X_{\mathcal{E}}+\lambda_{2}I_{|\mathcal{E}|})^{-1}$. Using Lemma 2.1, we get an algebraic curve $r_{j,\mathcal{E},s}(\lambda_{1},\lambda_{2})=0$ with degree 1 in $\lambda_{1}$ and $|\mathcal{E}|$ in $\lambda_{2}$ corresponding to addition of $j\notin\mathcal{E}$ given $\mathcal{E},s$.

Condition for leaving $\mathcal{E}$. Now consider a fixed $j^{\prime}\in\mathcal{E}$, given fixed $\mathcal{E},s$. The coefficient of $j^{\prime}$ will be zero for $\lambda_{1}^{*}=\frac{(c_{1})_{j^{\prime}}}{(c_{2})_{j^{\prime}}}$, which simplifies to $\lambda_{1}(({X_{\mathcal{E}}}^{T}X_{\mathcal{E}}+\lambda_{2}I_{|\mathcal{E}|})^{-1}s)_{j^{\prime}}=(({X_{\mathcal{E}}}^{T}X_{\mathcal{E}}+\lambda_{2}I_{|\mathcal{E}|})^{-1}{X_{\mathcal{E}}}^{T}y)_{j^{\prime}}$. Again by Lemma 2.1, we get an algebraic curve $r_{j^{\prime},\mathcal{E},s}(\lambda_{1},\lambda_{2})=0$ with degree 1 in $\lambda_{1}$ and at most $|\mathcal{E}|$ in $\lambda_{2}$ corresponding to removal of $j^{\prime}\in\mathcal{E}$ given $\mathcal{E},s$.

Putting the two together, we get $\sum_{i=0}^{p}2^{i}{p\choose i}\left((p-i)+i\right)=p3^{p}$ algebraic curves of degree 1 in $\lambda_{1}$ and at most $p$ in $\lambda_{2}$, across which the signed equicorrelation set may change. These curves characterize the complete set of points $(\lambda_{1},\lambda_{2})$ at which $(\mathcal{E},s)$ may possibly change. Thus by setting these $p3^{p}$ curves as the set of boundary functions $\mathcal{G}$, $\mathcal{E},s$ is guaranteed to be fixed for each sign pattern, and the corresponding loss takes the rational function form shown above. ∎

The exact same piecewise structure can be established for the dual function classes for loss functions $\ell_{\text{EN}}^{\textsf{AIC}}(\lambda,\cdot)$ and $\ell_{\text{EN}}^{\textsf{BIC}}(\lambda,\cdot)$. This is evident from the proof of Theorem 2.2, since any dual piece has a fixed equicorrelation set, and therefore $||\beta||_{0}$ is fixed. Given this piecewise structure, a challenge to learning values of $\lambda$ that minimize the loss function is that the function may not be differentiable (or may even be discontinuous, for the information criteria based losses) at the piece boundaries, making well-known gradient-based (local) optimization techniques inapplicable here. In the following (specifically Algorithm 1) we will show that techniques from data-driven design may be used to overcome this optimization challenge.

3 Learning to Regularize the ElasticNet

We will consider the problem of learning provably good ElasticNet parameters for a given problem domain, from multiple datasets (problem instances) either available as a collection (Section 3.1), or arriving online (Section 3.2). Our parameter tuning techniques also apply to simpler regression techniques typically used for variable selection, like LARS and LASSO, which are reasonable choices if the features are not multicollinear. Additional proof details for the results in this section are located in Appendix C.

3.1 Distributional Setting

Our main result in this section is the following upper bound on the pseudo-dimension of the classes of loss functions for the ElasticNet, which implies that in our distributional setting it is possible to learn near-optimal values of $\lambda$ with polynomially many problem instances.

Theorem 3.1.

$\textsc{Pdim}(\mathcal{H}_{\text{EN}})=O(p^{2})$. Further, $\textsc{Pdim}(\mathcal{H}_{\text{EN}}^{\textsf{AIC}})=O(p^{2})$ and $\textsc{Pdim}(\mathcal{H}_{\text{EN}}^{\textsf{BIC}})=O(p^{2})$.

Proof Sketch. The crucial ingredient is the $(\mathcal{F},\mathcal{G},p3^{p})$-piecewise decomposable structure for the dual class function $\mathcal{H}_{\text{EN}}^{*}$ established in Theorem 2.2, where $\mathcal{F}$ is a class of bivariate rational functions and $\mathcal{G}$ consists of polynomial thresholds, both with bounded degrees. We then bound the complexity of the corresponding dual class functions $\mathcal{F}^{*}$ and $\mathcal{G}^{*}$, in order to use the following powerful general result due to [BDD${}^{+}$21] (Theorem C.2 in the appendix)

$$\textsc{Pdim}(\mathcal{H})=O((\textsc{Pdim}(\mathcal{F}^{*})+d_{\mathcal{G}^{*}})\log(\textsc{Pdim}(\mathcal{F}^{*})+d_{\mathcal{G}^{*}})+d_{\mathcal{G}^{*}}\log k).$$

In more detail, we can bound the pseudo-dimension of the dual class of piece functions $\mathcal{F}^{*}$ (a class of bivariate rational functions) by $O(\log p)$ (Lemma C.4 in the appendix), by giving an upper bound of $O(k^{3}d^{3})$ on the number of sign patterns over $\mathbb{R}^{2}$ induced by $k$ algebraic curves of degree at most $d$. We can also bound the VC dimension of the dual class of boundary functions $\mathcal{G}^{*}$ (polynomial thresholds in two variates) by $O(p)$ using a standard linearization argument (Lemma C.5). Finally, the above result from [BDD${}^{+}$21] allows us to bound the pseudodimension of $\mathcal{H}$ by combining the above bounds.

$$\textsc{Pdim}(\mathcal{H})=O(p\log p+p\log(p3^{p}))=O(p^{2}).$$

The dual classes ${(\mathcal{H}_{\text{EN}}^{\textsf{AIC}})}^{*}$ and ${(\mathcal{H}_{\text{EN}}^{\textsf{BIC}})}^{*}$ also follow the same piecewise decomposable structure given by Theorem 2.2. This is because in each piece the equicorrelation set $\mathcal{E}$, and therefore $||\beta||_{0}=|\mathcal{E}|$ is fixed (Lemma B.2). The above argument implies an identical upper bound on the pseudo-dimensions of $\mathcal{H}_{\text{EN}}^{\textsf{AIC}}$ and $\mathcal{H}_{\text{EN}}^{\textsf{BIC}}$. See Appendix C for further proof details, including the technical lemmas. $\square$

The upper bound above implies a guarantee on the sample complexity of learning the ElasticNet tuning parameter, using standard learning-theoretic results [AB99], under mild boundedness assumptions on the data and hyperparameter search space.

Assumption 1 (Boundedness).

The predicted variable and all feature values are bounded by an absolute constant $R$, i.e. $\max\{||X^{(i)}||_{\infty,\infty},||y^{(i)}||_{\infty},||X_{\text{val}}^{(i)}||_{\infty,\infty},||y_{\text{val}}^{(i)}||_{\infty}\}\leq R$. Furthermore, the regularization coefficients are bounded, $(\lambda_{1},\lambda_{2})\in[\lambda_{\min},\lambda_{\max}]^{2}$ for $0<\lambda_{\min}<\lambda_{\max}<\infty$.

In our setting of learning from multiple problem instances, each sample is a dataset instance, so the sample complexity is simply the number of regression problem instances needed to learn the tuning parameters to any given approximation and confidence level.

Theorem 3.2 (Sample complexity of tuning the ElasticNet).

Suppose Assumption 1 holds. Let $\mathcal{D}$ be an arbitary distribution over the problem space $\Pi_{m,p}$. There is an algorithm which given $n=O\left(\frac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}H^{2}}}{\epsilon^{2}}(p^{2}+\log\frac{1}{\delta})\right)$ problem samples drawn from $\mathcal{D}$, for any $\epsilon>0$ and $\delta\in(0,1)$ and some constant $H$, outputs a regularization parameter $\hat{\lambda}$ for the ElasticNet such that with probability at least $1-\delta$ over the draw of the problem samples, we have that

$$\Big{\lvert}\mathbb{E}_{P\sim\mathcal{D}}[\ell_{EN}(\hat{\lambda},P)]-\min_{\lambda}\mathbb{E}_{P\sim\mathcal{D}}[\ell_{EN}(\lambda,P)]\Big{\rvert}\leq\epsilon.$$

Proof.

We use Lemma C.6 to conclude that the validation loss is uniformly bounded by some constant $H$ under Assumption 1. The result then follows from substituting our result in Theorem 3.1 into well-known generalization guarantee for function classes with bounded pseudo-dimensions (Theorem A.1). ∎

Discussion and applications. Computing the parameters which minimize the loss on the problem samples (aka Empirical Risk Minimization, or ERM) achieves the sample complexity bound in Theorem 3.2. Even though we only need polynomially many samples to guarantee the selection of nearly-optimal parameters, it is not clear how to implement the ERM efficiently. Note that we do not assume the set of features is the same across problem instances, so our approach can handle feature reset i.e. different problem instances can differ in not only the number of examples but also the number of features. Moreover, as a special case application, we consider the commonly used techniques of leave-one-out cross validation (LOOCV) and Monte Carlo cross validation (repeated random test-validation splits, typically independent and in a fixed proportion). Given a dataset of size $m_{tr}$, LOOCV would require $m_{tr}$ regression fits which can be inefficient for large dataset size. Alternately, we can consider draws from a distribution $\mathcal{D}_{LOO}$ which generates problem instances $P$ from a fixed dataset $(X,y)\in\mathbb{R}^{m+1\times p}\times\mathbb{R}^{m+1}$ by uniformly selecting $j\in[m+1]$ and setting $P=(X_{-j*},y_{-j},X_{j*},y_{j})$. Theorem 3.2 now implies that $\tilde{O}(p^{2}/\epsilon^{2})$ iterations are enough to determine an ElasticNet parameter $\hat{\lambda}$ with loss within $\epsilon$ (with high probability) of the parameter $\lambda^{*}$ obtained from running the full LOOCV. Similarly, we can define a distribution $\mathcal{D}_{MC}$ to capture the Monte Carlo cross validation procedure and determine the number of iterations sufficient to get an $\epsilon$-approximation of the loss corresponding parameter selection with arbitrarily large number of runs of the procedure. Thus, in a very precise sense, our results answer the question of how much cross-validation is enough to effectively implement the above techniques.

Remark 1.

While our result implies polynomial sample complexity, the question of learning the provably near-optimal parameter efficiently (even in output polynomial time) is left open. For the special cases of LASSO $(\lambda_{2}=0)$ and Ridge $(\lambda_{1}=0)$, the piece boundaries of the piecewise polynomial dual class (loss) function may be computed efficiently (using the LARS-LASSO algorithm of [EHJT04] for LASSO, and solving linear systems and locating roots of polynomials for Ridge). This applies to online and classification settings in the following sections as well.

3.2 Online Learning

We will now extend our results to learning the regularization coefficients given an online sequence of regression problems, such as when one needs to solve a new regression problem each day. Unlike the distributional setting above, we will not assume any problem distribution and our results will hold for an adversarial sequence of problem instances. We will need very mild assumptions on the data, namely boundedness of feature and prediction values and ‘smoothness’ of predictions (formally stated as Assumptions 1 and 2), while our distributional results above hold for worst-case problem datasets.

We will need two mild assumptions on the datasets in our problem instances for our results to hold. Our first assumption is that all feature values and predictions are bounded, for training as well as validation examples (Assumption 1 above). We will need the following definition to state our second assumption. Roughly speaking the definition below captures smoothness of a distribution.

Definition 3.

A continuous probability distribution is said to be $\kappa$-bounded if the probability density function $p(x)$ satisfies $p(x)\leq\kappa$ for any $x$ in the sample space.

For example, the normal distribution $\mathcal{N}(\mu,\sigma^{2})$ with mean $\mu$ and standard deviation $\sigma$ is $\frac{1}{\sigma\sqrt{2\pi}}$-bounded. We assume that the predicted variable $y$ in the training set comes from a $\kappa$-bounded (i.e. smooth) distribution, which does not require the strong tail decay of sub-Gaussian distributions [Zha09, CP09]. Moreover, the online adversary is allowed to change the distribution as long as it is $\kappa$-bounded. Note that our assumption also captures common data preprocessing steps, for example the jitter parameter in the popular Python library scikit-learn [P${}^{+}$11] adds a uniform noise to the $y$ values to help model stability. The assumption is formally stated as follows:

Assumption 2 (Smooth predictions).

The predicted variables $y^{(i)}$ in the training set are drawn from a joint $\kappa$-bounded distribution, i.e. for each $i$, the variables $y^{(i)}$ have a joint distribution with probability density bounded by $\kappa$.

Under these assumptions, we can show that it is possible to learn the ElasticNet parameters with sublinear expected regret when the problem instances arrive online. The learning algorithm (Algorithm 1) that achieves this regret is a continuous variant of the classic Exponential Weights algorithm [CBL06, BDV18]. It samples points in the domain with probability inversely propotional to the exponentiated loss. To formally state our result, we will need the following definition of dispersed loss functions. Informally speaking, it captures how amenable a set of non-Lipschitz functions is to online learning by measuring the worst rate of occurrence of non-Lipschitzness (or discontinuities) between any pair of points in the domain. [BDV18, BDS20, BDP20] show that dispersion is necessary and sufficient for learning piecewise Lipschitz functions.

Definition 4.

Dispersion [BDP20]. The sequence of random loss functions $l_{1},\dots,l_{T}$ is $\beta$-dispersed for the Lipschitz constant $L$ if, for all $T$ and for all $\epsilon\geq T^{-\beta}$, we have that, in expectation, at most $\tilde{O}(\epsilon T)$ functions (the soft-O notation suppresses dependence on quantities beside $\epsilon,T$ and $\beta$, as well as logarithmic terms) are not $L$-Lipschitz for any pair of points at distance $\epsilon$ in the domain $\mathcal{C}$. That is, for all $T$ and for all $\epsilon\geq T^{-\beta}$, $\mathbb{E}\Big{[}\max_{\begin{subarray}{c}\rho,\rho^{\prime}\in\mathcal{C}\\ \left\lVert\rho-\rho^{\prime}\right\rVert_{2}\leq\epsilon\end{subarray}}\big{\lvert}\{t\in[T]\mid l_{t}(\rho)-l_{t}(\rho^{\prime})>L\left\lVert\rho-\rho^{\prime}\right\rVert_{2}\}\big{\rvert}\Big{]}\leq\tilde{O}(\epsilon T)$.

Our key contribution is to show that the loss sequence is dispersed (Definition 4) under the above assumptions. This involves establishing additional structure for the problem, specifically about the location of boundary functions in the piecewise structure from Theorem 2.2. This stronger characterization coupled with results from [BDP20] on dispersion of algebraic discontinuities completes the proof.

Theorem 3.3.

Suppose Assumptions 1 and 2 hold. Let $l_{1},\dots,l_{T}:(0,\lambda_{\max})^{2}\rightarrow\mathbb{R}_{\geq 0}$ denote an independent sequence of losses (e.g. fresh randomness is used to generate the validation set features in each round) as a function of the ElasticNet regularization parameter $\lambda=(\lambda_{1},\lambda_{2})$, $l_{i}(\lambda)=l_{r}(\hat{\beta}^{(X^{(i)},y^{(i)})}_{\lambda,f_{EN}},(X_{\text{val}}^{(i)},y_{\text{val}}^{(i)}))$. The sequence of functions is $\frac{1}{2}$-dispersed, and there is an online algorithm with $\tilde{O}(\sqrt{T})$ expected regret. The result also holds for loss functions adjusted by information criteria AIC and BIC.

Proof Sketch. We start with the $(\mathcal{F},\mathcal{G},p3^{p})$-piecewise decomposable structure for the dual class function $\mathcal{H}_{\text{EN}}^{*}$ from Theorem 2.2. Observe that the rational piece functions in $\mathcal{F}$ do not introduce any new discontinuities since the denominator polynomials do not have positive roots. For each of two types of boundary functions in $\mathcal{G}$ (corresponding to leaving/entering the equicorrelation set) we show that the discontinuities between any pair of points $\lambda,\lambda^{\prime}$ lie along the roots of polynomials with non-leading coefficients bounded and smoothly distributed (bounded joint density). This allows us to use results from [BDP20] to establish dispersion, and therefore online learnability. $\square$

1:  Input: Problems $(X^{(i)},y^{(i)})$ and regularization penalty function $f$.

2:  Hyperparameter: step size parameter $\zeta\in(0,1]$.

3:  Output: Regularization parameter $(\lambda_{i})_{i\in[T]}\in C$, $C\subset\mathbb{R}^{+}$ (LASSO/Ridge) or $C\subset{\mathbb{R}^{+}}^{2}$ (ElasticNet).

4:  Set $w_{1}(\lambda)=1$ for all $\lambda\in C$.

5:  for $i=1,2,\dots,T$ do

6:     $W_{i}:=\int_{C}w_{i}(\lambda)d\lambda$.

7:     Sample $\lambda$ with probability $p_{t}(\lambda)=\frac{w_{i}(\lambda)}{W_{i}}$, output as $\lambda_{i}$.

8:     Compute average loss function $l_{i}(\lambda)=\frac{1}{|y^{(i)}|}l(\hat{\beta}_{\lambda,f},(X^{(i)},y^{(i)}))$.

9:     For each $\lambda\in C,\text{ update weights }w_{i+1}(\lambda)=e^{\zeta(1-l_{i}(\lambda))}w_{i}(\lambda)$.

10:  end for

Algorithm 1 Data-driven Regularization ($\zeta$)

We remark that the above result holds for arbitrary training features and validation sets in the problem sequence that satisfy our assumptions, in particular the losses are only assumed to be independent but not identically distributed. In contrast, the results in the previous section needed them to be drawn from the same distribution. Also the parameters need to be selected online, and cannot be changed for already seen instances. This setting captures interesting practical settings where the set of features (including feature dimensions) and the relevant training set (including training set size) may change over the online sequence. It is not clear how usual model selection techniques like cross-validation may be adapted to these challenging settings.

4 Extension to Regularized Least Squares Classification

Regression techniques can also be used to train binary classifiers by using an appropriate threshold on top of the regression estimate. Intuitively, regression learns a linear mapping which projects the datapoints onto a one-dimensional space, i.e. a real number, after which a threshold may be applied to classify the points. The use of thresholds to make discrete classifications adds discontinuities to the empirical loss function. Thus, in general, the classification setting is more challenging as it already includes the piecewise structure in the regression loss. We provide statistical and online learning guarantees for Ridge and LASSO. For the ElasticNet we present the extensions needed to the arguments from the previous sections to obtain results in the classification setting.

More formally, we will restrict $y$ to $\{0,1\}^{m}$. The estimator $\hat{\beta}_{\lambda,f}$ is obtained as before, and the prediction on a test instance $x$ may be obtained by taking the sign of a thresholded regression estimate, $\textup{{sgn}}(\langle x,\hat{\beta}_{\lambda,f}\rangle-\tau)$, where $\textup{{sgn}}:\mathbb{R}\rightarrow\{0,1\}$ maps $x\in\mathbb{R}$ to $\mathbb{I}\{x\geq 0\}$ and $\tau\in\mathbb{R}$ is the threshold. The threshold $\tau$ corresponds to the intercept or bias of the learned linear classifier, here we will treat it as a tunable hyperparameter (in addition to $\lambda_{1},\lambda_{2}$). The average 0-1 loss over the dataset $(X,y)$ is given by $l_{c}(\hat{\beta}_{\lambda,f},(X,y),\tau)=\frac{1}{m}\sum_{i=1}^{m}|y_{i}-\textup{{sgn}}(\langle X_{i},\hat{\beta}_{\lambda,f}\rangle-\tau)|$. Proofs from this section are in Appendix D.

4.1 Distributional setting

The problem setting is the same as in Section 3.1, except that the labels $y$ are binary and we use threshold for prediction. We bound the pseudo-dimension for classification loss on these problem instances, which as before (c.f. Theorems 3.1 and 3.2) imply that polynomially many problem samples are sufficient to generalize well over the problem distribution $\mathcal{D}$. For Ridge and LASSO we upper bound the number of discontinuities of the piecewise constant classification loss by determining the values of $\lambda$ where any prediction changes.

Theorem 4.1.

Let $\mathcal{H}_{\text{Ridge}}^{c}$, $\mathcal{H}_{\text{LASSO}}^{c}$ and $\mathcal{H}_{\text{EN}}^{c}$ denote the set of loss functions for classification problems with at most $m$ examples and $p$ features, for linear classifiers regularized using Ridge, LASSO and ElasticNet regression respectively.

• $\textsc{Pdim}(\mathcal{H}_{\text{Ridge}}^{c})=O(\log mp)$

• $\textsc{Pdim}(\mathcal{H}_{\text{LASSO}}^{c})=O(p\log m)$. Further, in the overparameterized regime ($p\gg m$), we have that $\textsc{Pdim}(\mathcal{H}_{\text{LASSO}}^{c})=O(m\log\frac{p}{m})$.

• $\textsc{Pdim}(\mathcal{H}_{\text{EN}}^{c})=O(p^{2}+p\log m)$.

The key difference with the bound for the regression loss in Theorem 3.1 is the additional $O(p\log m)$ term which corresponds to discontinuities induced by the thresholding in the regression based classifiers. We can establish a structure similar to Theorem 2.2 in this case (Lemma D.1).

4.2 Online setting

As in Section 3.2, we can define an online learning setting for classification. Note that the smoothness of the predicted variable is not meaningful here, since $y$ is a binary vector. Instead we will assume that the validation examples have smooth feature values. Intuitively this means that small perturbations to the feature values does not meaningfully change the problem.

Assumption 3 (Smooth validation features).

The feature values $(X_{\text{val}}^{(i)})_{jk}$ in the validation examples are drawn from a joint $\kappa$-bounded distribution.

Under the assumption, we show that we can learn the regularization parameters online, for each of Ridge, LASSO and ElasticNet estimators. The proofs are straightforward extensions of the structural results developed in the previous sections, with minor technical changes to use the above validation set feature smoothness instead of Assumption 2, and are deferred to the appendix.

Theorem 4.2.

Suppose Assumptions 1 and 3 hold. Let $l_{1},\dots,l_{T}:(0,H]^{d}\times[-H,H]\rightarrow\mathbb{R}$ denote an independent sequence of losses as a function of the regularization parameter $\lambda$, $l_{i}(\lambda,\tau)=l_{c}(\hat{\beta}_{\lambda,f},(X^{(i)},y^{(i)}),\tau)$. If $f$ is given by $f_{1}$ (LASSO), $f_{2}$ (Ridge), or $f_{EN}$ (ElasticNet) then the sequence of functions is $\frac{1}{2}$-dispersed and there is an online algorithm with $\tilde{O}(\sqrt{T})$ expected regret.

5 Conclusions and Future Work

We obtain a novel structural result for the ElasticNet loss as a function of the tuning parameters. Our characterization is useful in giving upper bounds for the sample complexity of learning the parameters from multiple regression problem instances (i.e. different datasets, possibly corresponding to different tasks) from the same problem domain. Efficient algorithms are immediate from our results for Ridge and LASSO. For the ElasticNet we show generalization and online regret guarantees, but efficient implementation of the algorithms is an interesting question for further work. Also we show general learning-theoretic guarantees, i.e. without any significant restrictions on the data-generating distribution, in learning from multiple problems. The problems may be drawn i.i.d. from an arbitrary problem distribution, or even arrive in an online sequence but with some smoothness properties. It is unclear if such general guarantees may be given for tuning parameters for the more standard setting of tuning over a single training set generated by i.i.d. draws from an example distribution, or how such guarantees can be combined with our results.

Appendix

Appendix A A classic Generalization Bound

The pseudo-dimension (also known as the Pollard dimension) is a generalization of the VC-dimension to real-valued functions, and may be defined as follows.

Definition 5 (Pseudo-dimension [Pol12]).

Let $\mathcal{H}$ be a set of real valued functions from input space $\mathcal{X}$. We say that $C=(x_{1},\dots,x_{n})\in\mathcal{X}^{n}$ is pseudo-shattered by $\mathcal{H}$ if there exists a vector $r=(r_{1},\dots,r_{n})\in\mathbb{R}^{n}$ (called “witness”) such that for all $b=(b_{1},\dots,b_{n})\in\{\pm 1\}^{n}$ there exists $h_{b}\in\mathcal{H}$ such that $\text{sign}(h_{b}(x_{i})-r_{i})=b_{i}$. Pseudo-dimension of $\mathcal{H}$, denoted by $\textsc{Pdim}(\mathcal{H})$, is the cardinality of the largest set pseudo-shattered by $\mathcal{H}$.

The following theorem connects the sample complexity of uniform learning over a class of real-valued functions to the pseudo-dimension of the class. Let $h^{*}:\mathcal{X}\rightarrow\{0,1\}$ denote the target concept. We say $\mathcal{H}$ is $(\epsilon,\delta)$-uniformly learnable with sample complexity $n$ if, for every distribution $\mathcal{D}$, given a sample $S$ of size $n$, with probability $1-\delta$, $\big{\lvert}\frac{1}{n}\sum_{s\in S}|h(s)-h^{*}(s)|-\mathbb{E}_{s\sim\mathcal{D}}[|h(s)-h^{*}(s)|]\big{\rvert}<\epsilon$ for every $h\in\mathcal{H}$.

Theorem A.1 ([AB99]).

Suppose $\mathcal{H}$ is a class of real-valued functions with range in $[0,H]$ and pseudo-dimension $\textsc{Pdim}(\mathcal{H})$. For every $\epsilon>0,\delta\in(0,1)$, the sample complexity of $(\epsilon,\delta)$-uniformly learning the class $\mathcal{H}$ is $O\left(\left(\frac{H}{\epsilon}\right)^{2}\left(\textsc{Pdim}(\mathcal{H})\ln\frac{H}{\epsilon}+\ln\frac{1}{\delta}\right)\right)$.

Appendix B Known characterization of LASSO solutions

We will review some properties of LASSO solutions from prior work that are useful in proving our results. Let $(X,y)$ with $X=[\bm{x}_{1},\dots,\bm{x}_{p}]\in\mathbb{R}^{m\times p}$ and $y\in\mathbb{R}^{m}$ denote a (training) dataset consisting of $m$ labeled examples with $p$ features. As noted in Section 2, LASSO regularization may be formulated as the following optimization problem.

$$\min_{\beta\in\mathbb{R}^{p}}\left\lVert y-X\beta\right\rVert_{2}^{2}+\lambda_{1}||\beta||_{1},$$

where $\lambda_{1}\in\mathbb{R}^{+}$ is the regularization parameter. Dealing with the case $\lambda_{1}=0$ (i.e. Ordinary Least Squares) is not difficult, but is omitted here to keep the statements of the definitions and results simple. We will use the following well-known facts about the solution of the LASSO optimization problem [Fuc05, Tib13]. Applying the Karush-Kuhn-Tucker (KKT) optimality conditions to the problem gives,

Lemma B.1 (KKT Optimality Conditions for LASSO).

$\beta^{*}\in\operatorname*{argmin}_{\beta\in\mathbb{R}^{p}}\left\lVert y-X\beta\right\rVert_{2}^{2}+\lambda_{1}||\beta||_{1}$ iff for all $j\in[p]$,

$$\bm{x}_{j}^{T}(y-X\beta^{*}) =\lambda_{1}\textup{{sign}}(\beta^{*}),\text{ if }\beta^{*}_{j}\neq 0,$$

$$|\bm{x}_{j}^{T}(y-X\beta^{*})| \leq\lambda_{1},\text{ otherwise. }$$

Here $\bm{x}_{j}^{T}(y-X\beta^{*})$ is simply the correlation of the the $j$-th covariate with the residual $y-X\beta^{*}$ (when $y,X$ have been standardized). This motivates the definition of equicorrelation sets of covariates (Definition 2).

In terms of the equicorrelation set and the equicorrelation sign vector, the characterization of the LASSO solution in Lemma B.1 implies

$$X_{\mathcal{E}}^{T}(y-X_{\mathcal{E}}\beta^{*}_{\mathcal{E}})=\lambda_{1}s.$$

This implies a necessary and sufficient condition for the uniqueness of the LASSO solution, namely that $X_{\mathcal{E}}$ is full rank for all equicorrelation sets $\mathcal{E}$ [Tib13]. Our results will hold if the dataset $X$ satisfies this condition, but for simplicity we will use the a simpler (and possibly more natural) sufficient condition involving the general position.

Definition 6.

A matrix $X\in\mathbb{R}^{m\times p}$ is said to have its columns in the general position if the affine span of any $k\leq m$ points $(\sigma_{i}\bm{x}_{j_{i}})_{i\in[k],\{j_{i}\}_{i}=J\subseteq[p]}$ for arbitrary signs $\sigma_{[k]}\in\{-1,1\}^{k}$ and subset $J$ of the columns of size $k$, does not contain any element of $\{\bm{x}_{i}\mid i\notin J\}$.

Finally, we state the following useful characterization of the LASSO solutions in terms of the equicorrelation sets and sign vectors.

Lemma B.2 ([Tib13], Lemma 3).

If the columns of $X$ are in general position, then for any $y$ and $\lambda_{1}>0$, the LASSO solution is unique and is given by

$$\beta^{*}_{\mathcal{E}}=(X_{\mathcal{E}}^{T}X_{\mathcal{E}})^{-1}(X_{\mathcal{E}}^{T}y-\lambda_{1}s),\beta^{*}_{[p]\setminus\mathcal{E}}=0.$$

We remark that Lemma B.2 does not give a way to compute $\beta^{*}$ for a given value of $\lambda_{1}$, since $\mathcal{E}$ and $s$ depend on $\beta^{*}$, but still gives a property of $\beta^{*}$ that is convenient to use. In particular, since we have at most $3^{p}$ possible choices for $(\mathcal{E},s)$, this implies that the LASSO solution $\beta^{*}(\lambda_{1})$ is a piecewise linear function of $\lambda_{1}$, with at most $3^{p}$ pieces (for $\lambda_{1}>0$). Following popular terminology, we will refer to this function as a solution path of LASSO for the given dataset $(X,y)$. LARS-LASSO of [ZH05] is an efficient algorithm for computing the solution path of LASSO.

Corollary B.3.

Let $X$ be a matrix with columns in the general position. If the unique LASSO solution for the dataset $(X,y)$ is given by the function $\beta^{*}:\mathbb{R}^{+}\rightarrow\mathbb{R}^{p}$, then $\beta^{*}$ is piecewise linear with at most $3^{p}$ pieces given by Lemma B.2.

Appendix C Lemmas and proof details for Section 3

We start with a helper lemma that characterizes the solution of the ElasticNet in terms of equicorrelation sets and sign vectors.

Lemma C.1.

Let $X$ be a matrix with columns in the general position, and $\lambda=(\lambda_{1},\lambda_{2})\in(0,\infty)\times(0,\infty)$. Then the ElasticNet solution $\hat{\beta}_{\lambda,f_{EN}}\in\operatorname*{argmin}_{\beta\in\mathbb{R}^{p}}\left\lVert y-X\beta\right\rVert_{2}^{2}+\langle\lambda,f_{\text{EN}}(\beta)\rangle$ is unique for any dataset $(X,y)$ and satisfies

$$\hat{\beta}_{\lambda,f_{EN}}=(X_{\mathcal{E}}^{T}X_{\mathcal{E}}+\lambda_{2}I_{|\mathcal{E}|})^{-1}X_{\mathcal{E}}^{T}y-\lambda_{1}(X_{\mathcal{E}}^{T}X_{\mathcal{E}}+\lambda_{2}I_{|\mathcal{E}|})^{-1}s$$

for some $\mathcal{E}\subseteq[p]$ and $s\in\{-1,1\}^{|\mathcal{E}|}$.

Proof.

We start with the well-known characterization of the ElasticNet solution as the solution of a LASSO problem on a transformed dataset, obtained using simple algebra [ZH05]. Given any dataset $(X,y)$, the ElasticNet coefficients $\hat{\beta}_{\lambda,f_{EN}}$ are given by $\hat{\beta}_{\lambda,f_{EN}}=\frac{1}{\sqrt{1+\lambda_{2}}}\hat{\beta}^{*}_{\lambda}$ where $\hat{\beta}^{*}_{\lambda}$ is the solution for a LASSO problem on a modified dataset $(X^{*},y^{*})$

$$\hat{\beta}^{*}_{\lambda}=\operatorname*{argmin}_{\beta}\left\lVert y^{*}-X^{*}\beta\right\rVert_{2}^{2}+\lambda_{1}^{*}f_{1}(\beta),$$

with $X^{*}=\frac{1}{\sqrt{1+\lambda_{2}}}\begin{pmatrix}X\\ \sqrt{\lambda_{2}}I_{p}\end{pmatrix}$, $y^{*}=\begin{pmatrix}y\\ 0\end{pmatrix}$, and $\lambda_{1}^{*}=\frac{\lambda_{1}}{\sqrt{1+\lambda_{2}}}$.

If the columns of $X$ are in general position (Definition 6), then the same is true of $X^{*}$. For $\mathcal{E}\subseteq[p]$, note that ${X^{*}_{\mathcal{E}}}^{T}X^{*}_{\mathcal{E}}=\frac{1}{1+\lambda_{2}}(X_{\mathcal{E}}^{T}X_{\mathcal{E}}+\lambda_{2}I_{|\mathcal{E}|})$ and ${X^{*}_{\mathcal{E}}}^{T}y^{*}=\frac{1}{\sqrt{1+\lambda_{2}}}X_{\mathcal{E}}^{T}y$. By Lemma B.2, if $\mathcal{E}$ denotes the equicorrelation set of covariates and $s\in\{-1,1\}^{|\mathcal{E}|}$ the equicorrelation sign vector for the LASSO problem, then the ElasticNet solution is given by

$$\hat{\beta}_{\lambda,f_{EN}}=c_{1}-c_{2}{\lambda_{1}},$$

where

$$c_{1}=\frac{1}{\sqrt{1+\lambda_{2}}}({X^{*}_{\mathcal{E}}}^{T}X^{*}_{\mathcal{E}})^{-1}{X^{*}_{\mathcal{E}}}^{T}y^{*}=(X_{\mathcal{E}}^{T}X_{\mathcal{E}}+\lambda_{2}I_{|\mathcal{E}|})^{-1}X_{\mathcal{E}}^{T}y,$$

and

$$c_{2}=\frac{1}{1+\lambda_{2}}({X^{*}_{\mathcal{E}}}^{T}X^{*}_{\mathcal{E}})^{-1}s=(X_{\mathcal{E}}^{T}X_{\mathcal{E}}+\lambda_{2}I_{|\mathcal{E}|})^{-1}s.$$

∎

The following lemma helps determine the dependence of ElasticNet solutions on $\lambda_{2}$.

Lemma 2.1 (restated). Let $A$ be an $r\times s$ matrix. Consider the matrix $B(\lambda)=(A^{T}A+\lambda I_{s})^{-1}$ for $\lambda>0$.

• Each entry of $B(\lambda)$ is a rational polynomial $P_{ij}(\lambda)/Q(\lambda)$ for $i,j\in[s]$ with each $P_{ij}$ of degree at most $s-1$, and $Q$ of degree $s$.

• Further, for $i=j$, $P_{ij}$ has degree $s-1$ and leading coefficient 1, and for $i\neq j$ $P_{ij}$ has degree at most $s-2$. Also, $Q(\lambda)$ has leading coefficient $1$.

Proof.

Let $G=A^{T}A$ be the Gramian matrix. $G$ is symmetric and therefore diagonalizable, and the diagonalization gives the eigendecomposition $G=E\Lambda E^{-1}$. Thus we have

$$(A^{T}A+\lambda I_{s})^{-1}=(E\Lambda E^{-1}+\lambda EE^{-1})^{-1}=E(\Lambda+\lambda I_{s})^{-1}E^{-1}$$

But $\Lambda$ is the diagonal matrix $\textsf{Diag}(\Lambda_{11},\dots,\Lambda_{ss})$, and therefore $(\Lambda+\lambda I_{s})^{-1}=\textsf{Diag}((\Lambda_{11}+\lambda)^{-1},\dots,(\Lambda_{ss}+\lambda)^{-1})$. This implies the desired characterization, with $Q(\lambda)=\Pi_{i\in[s]}(\Lambda_{ii}+\lambda)$ and

$$P_{ij}(\lambda)=Q(\lambda)\sum_{k=1}^{s}\frac{E_{ik}(E^{-1})_{kj}}{\Lambda_{kk}+\lambda}=\sum_{k=1}^{s}\left(E_{ik}(E^{-1})_{kj}\Pi_{i\in[s]\setminus k}(\Lambda_{ii}+\lambda)\right),$$

with coefficient of $\lambda^{s-1}$ in $P_{ij}(\lambda)$ equal to $\sum_{k=1}^{s}E_{ik}(E^{-1})_{kj}=\mathbb{I}\{i=j\}$. ∎

C.1 Tuning the ElasticNet – Distributional setting

We first present some terminology from algebraic geometry which will be useful in our proofs.

Definition 7 (Semialgebraic sets, Algebraic curves.).

A semialgebraic subset of $\mathbb{R}^{n}$ is a finite union of sets of the form $\{x\in\mathbb{R}^{n}\mid p_{i}(x)\geq 0\text{ for each }i\in[m]\}$, where $p_{1},\dots,p_{m}$ are polynomials. An algebraic curve is the zero set of a polynomial in two dimensions.

The result of Theorem 2.2 motivates the following results for bounding the complexity of dual piece functions and dual boundary functions, which can be used to bound the pseudo-dimension of $\mathcal{H}_{\text{EN}}$ (Theorem 3.1) using the following remarkable result from [BDD${}^{+}$21].

Theorem C.2 ([BDD${}^{+}$21]).

If the dual function class $\mathcal{H}^{*}$ is $(\mathcal{F},\mathcal{G},k)$-piecewise decomposable, then the pseudo-dimension of $\mathcal{H}$ may be bounded as

$$\textsc{Pdim}(\mathcal{H})=O((\textsc{Pdim}(\mathcal{F}^{*})+d_{\mathcal{G}^{*}})\log(\textsc{Pdim}(\mathcal{F}^{*})+d_{\mathcal{G}^{*}})+d_{\mathcal{G}^{*}}\log k),$$

where $d_{\mathcal{G}^{*}}$ denotes the VC dimension of dual class boundary function $\mathcal{G}^{*}$.

We will first prove a useful lemma that bounds the number of pieces into which a finite set of algebraic curves with bounded degrees may partition $\mathbb{R}^{2}$.

Lemma C.3.

Let $\mathcal{H}$ be a collection of $k$ functions $h_{i}:\mathbb{R}^{2}\rightarrow\mathbb{R}$ that map $(x,y)\mapsto q_{i}(x,y)$ where $q_{i}$ is a bivariate polynomial of degree at most $d$, for $i\in[k]$, then $\mathbb{R}^{2}\setminus\{(x,y)\mid q_{i}(x,y)=0\text{ for some }i\in[k]\}$ may be partitioned into at most $(kd+1)\left(d^{2}{k\choose 2}+2kd(d-1)+1\right)=O(d^{3}k^{3})$ disjoint sets such that the sign pattern $(\mathbb{I}\{q_{i}(x,y)>0\})_{i\in[k]}$ is fixed over any set in the partition.

Proof.

Assume WLOG that the curves are in the general position. Simple applications of Bezout’s theorem (which states that, in general, two algebraic curves of degrees $m$ and $n$ intersect in at most $mn$ points) imply that there are at most $d^{2}{k\choose 2}$ points where any pair of curves from the set $\{q_{i}(x,y)\}_{i\in[k]}$ may intersect, and at most $2kd(d-1)$ points of extrema (i.e. points $p_{0}=(x_{0},y_{0})$ on the curve $f$ such that there is an open neighborhood $N$ around $p_{0}$ in which $x_{0}\in\operatorname*{argmin}_{(x,y)\in N\cap f}x$, or $x_{0}\in\operatorname*{argmax}_{(x,y)\in N\cap f}x$, or $y_{0}\in\operatorname*{argmin}_{(x,y)\in N\cap f}y$, or $y_{0}\in\operatorname*{argmax}_{(x,y)\in N\cap f}y$) for the $k$ algebraic curves. Let $\mathcal{P}$ denote the set of these $\leq d^{2}{k\choose 2}+2kd(d-1)$ points.

Now a horizontal line $y=c$ will have the exact same set of intersections with all the curves in $\mathcal{H}$ as a line $y=c^{\prime}$, and in the same order (including multiplicities), if none of the points in $\mathcal{P}$ lie between these lines. There are thus at most $|\mathcal{P}|+1$ distinct sequences of the $k$ curves that may correspond to the intersection sequence of any horizontal line. Moreover, any such horizontal line may intersect any curve in the set at most $d$ times (since a polynomial in degree $d$ has at most $d$ zeros), or at most $kd$ intersections with all the curves. Summing up over the distinct intersection sequences, we have at most $(kd+1)(|\mathcal{P}|+1)$ distinct sign patterns induced by the set of curves. ∎

We will now use Lemma C.3 to bound the pseudo-dimension of the relevant function classes (Theorem 2.2).

Lemma C.4.

Let $\mathcal{F}^{*}=\{f^{*}_{q_{1},q_{2}}:\mathbb{R}^{2}\rightarrow\mathbb{R}\}$ be a function class consisting of rational polynomial functions $f^{*}_{q_{1},q_{2}}:(\lambda_{1},\lambda_{2})\mapsto\frac{q_{1}(\lambda_{1},\lambda_{2})}{q_{2}(\lambda_{1},\lambda_{2})}$, where $q_{1},q_{2}$ have degrees at most $d$. Then $\textsc{Pdim}(\mathcal{F}^{*})=O(\log d)$.

Proof.

Suppose that $\textsc{Pdim}(\mathcal{F}^{*})=N$. Then there exist functions $f^{*}_{1},\dots,f^{*}_{N}\in\mathcal{F}^{*}$ and real-valued witnesses $(r_{1},\dots,r_{N})\in\mathbb{R}^{N}$ such that for every subset $T\subseteq[N]$, there exists a parameter setting $\lambda_{T}=(\lambda_{1},\lambda_{2})\in\mathbb{R}^{2}$ such that $f^{*}_{i}(\lambda_{T})\geq r_{i}$ if and only if $i\in T$. In other words, we have a set of $2^{N}$ parameters (indexed by $T$) that induce all possible labelings of the binary vector $(\mathbb{I}\{f^{*}_{i}(\lambda_{T})\geq r_{i}\})_{i\in[N]}$.

But $f^{*}_{i}(\lambda)\geq r_{i}$ are semi-algebraic sets bounded by $N$ algebraic curves of degree at most $d$. By Lemma C.3, there are at most $O(d^{3}N^{3})$ different sign-patterns induced by $N$ algebraic curves over all possible values of $\lambda\in\mathbb{R}^{2}$. In particular, the number of distinct sign patterns over $\lambda\in\{\lambda_{T}\}_{T\subseteq[N]}$ is also $O(d^{3}N^{3})$. Thus, we conclude $2^{N}=O(d^{3}N^{3})$, or $N=O(\log d)$. ∎

Lemma C.5.

Let $R[x_{1},x_{2},\dots,x_{d}]_{D}$ denote the set of all real polynomials in $d$ variables of degree at most $D$ in $x_{1}$, and degree at most 1 in $x_{2},\dots,x_{d}$. Further, let $P_{d,D}=\{\{x\in R^{d}:p(x)\geq 0\}\mid p\in R[x_{1},x_{2},\dots,x_{d}]_{D}\}$. The VC-dimension of the set system $(\mathbb{R}^{d},P_{d,D})$ is $O(dD)$.

Proof.

We will employ a standard linearization argument [Cov65] that reduces the problem to bounding the VC dimension of halfspaces in higher dimensions. Let $M$ be the set of all possible non-constant monomials of degree at most $D$ in $x_{1}$, and at most one in $x_{2},\dots,x_{d}$. For example, when $d=3$ and $D=2$, we have $M=\{x_{1},x_{2},x_{3},x_{1}x_{2},x_{1}x_{3},x_{1}^{2},x_{1}^{2}x_{2},x_{1}^{2}x_{3}\}$. Note that $|M|=(D+1)d-1$. Indeed for each $x_{1}^{i}$ for $i=0,\dots,D$ we obtain a monomial by multiplying with each of $\{1,x_{2},\dots,x_{d}\}$. Excluding the constant monomial gives the result. The linearization we use is a map $\phi:\mathbb{R}^{d}\rightarrow\mathbb{R}^{|M|}$ which indexes the coordinates by monomials in $M$. For example when $d=3$ and $D=2$, $\phi(x_{1},x_{2},x_{3})=(x_{1},x_{2},x_{3},x_{1}x_{2},x_{1}x_{3},x_{1}^{2},x_{1}^{2}x_{2},x_{1}^{2}x_{3})$.

Now, if $S\in\mathbb{R}^{d}$ is shattered by $P_{d,D}$, then $\phi(S)$ is shattered by half-spaces in $\mathbb{R}^{|M|}$. Indeed, suppose $p=p_{0}+\langle\mathbf{p},\phi(x_{1},\dots,x_{d})\rangle\in P_{d,D}$ (for $\mathbf{p}\in\mathbb{R}^{|M|}$) is a polynomial that is positive over some $T\subseteq S$ and negative over $S\setminus T$. Define halfspace $h_{p}\in\mathbb{R}^{|M|}$ as $\{y\in\mathbb{R}^{|M|}\mid p_{0}+\langle\mathbf{p},y\rangle\geq 0\}$. Clearly $h_{p}\cap\phi(S)=\phi(T)$, and in general $S$ is shattered by halfspaces in $\mathbb{R}^{|M|}$. Using the well-known result for the VC-dimension of halfspaces we have that the VC-dimension of $P_{d,D}$ over $\mathbb{R}^{d}$ is $(D+1)d$. ∎

Theorem 3.1 (restated). $\textsc{Pdim}(\mathcal{H}_{\text{EN}})=O(p^{2})$. Further, $\textsc{Pdim}(\mathcal{H}_{\text{EN}}^{\textsf{AIC}})=O(p^{2})$ and $\textsc{Pdim}(\mathcal{H}_{\text{EN}}^{\textsf{BIC}})=O(p^{2})$.

Proof.

By Theorem 2.2, the dual class $\mathcal{H}_{\text{EN}}^{*}$ of $\mathcal{H}_{\text{EN}}$ is $(\mathcal{F},\mathcal{G},p3^{p})$-piecewise decomposable, with $\mathcal{F}=\{f_{q_{1},q_{2}}:\mathcal{L}\rightarrow\mathbb{R}\}$ consisting of rational polynomial functions $f_{q_{1},q_{2}}:l_{\lambda}\mapsto\frac{q_{1}(\lambda_{1},\lambda_{2})}{q_{2}(\lambda_{2})}$, where $q_{1},q_{2}$ have degrees at most $2p$, and $\mathcal{G}=\{g_{r}:\mathcal{L}\rightarrow\{0,1\}\}$ consisting of semi-algebraic sets bounded by algebraic curves $g_{r}:u_{\lambda}\mapsto\mathbb{I}\{r(\lambda_{1},\lambda_{2})<0\}$, where $r$ is a polynomial of degree 1 in $\lambda_{1}$ and at most $p$ in $\lambda_{2}$.

Now by Lemma C.4, we have $\textsc{Pdim}(\mathcal{F}^{*})=O(\log p)$, and by Lemma C.5 the VC dimension of the dual boundary class is $d_{\mathcal{G}^{*}}=O(p)$. A straightforward application of Theorem C.2 yields

$$\textsc{Pdim}(\mathcal{H})=O(p\log p+p\log(p3^{p}))=O(p^{2}).$$

The dual classes ${(\mathcal{H}_{\text{EN}}^{\textsf{AIC}})}^{*}$ and ${(\mathcal{H}_{\text{EN}}^{\textsf{BIC}})}^{*}$ also follow the same piecewise decomposable structure given by Theorem 2.2. This is because in each piece the equicorrelation set $\mathcal{E}$, and therefore $||\beta||_{0}=|\mathcal{E}|$ (by Lemma B.2) is fixed. Thus we can keep the same boundary functions, and the function value in each piece only changes by a constant (in $\lambda$) and is therefore also a rational function with the same degrees. The above argument then implies an identical upper bound on the pseudo-dimensions. ∎

The following lemma shows that under mild boundedness assumptions on the data and the search space of hyperparameters, the validation loss function class $\mathcal{H}_{\text{EN}}$ is uniformly bounded by some constant $H>0$.

Lemma C.6.

Under Assumption 1, there exists a constant $H>0$ so that for all $h(\lambda,\cdot)\in\mathcal{H}_{\text{EN}}=\{h(\lambda,\cdot):\Pi_{m,p}\rightarrow\mathbb{R}_{\geq 0}\mid\lambda\in[\lambda_{\min},\lambda_{\max}]\}$, we have $\|h(\lambda,\cdot)\|_{\infty}=\sup_{\begin{subarray}{c}P\in\Pi_{m,p}\end{subarray}}{h(\lambda,P)}\leq H$.

Proof.

For any problem instance $P=(X,y,X_{\text{val}},y_{\text{val}})\in\Pi_{m,p}$, and for any $\lambda=(\lambda_{1},\lambda_{2})\in[\lambda_{\min},\lambda_{\max}]^{2}$, consider the optimization problem for training set $(X,y)$

$$\operatorname*{argmin}_{\beta}F(\beta),$$ (2)

where $F(\beta)=\frac{1}{2m}||{y-X\beta}||_{2}^{2}+\lambda_{1}||\beta||_{1}+\lambda_{2}||\beta||_{2}^{2}$. If we set $\beta=\vec{0}$, we have

$$F(\vec{0})=\frac{1}{2m}||{y}||_{2}^{2}\leq\frac{1}{2}R^{2},$$

for some absolute constant $R$, due to Assumption 1. Let $\hat{\beta}_{(X,y)}(\lambda)$ be the optimal solution of 2, we have

$$\frac{1}{2}R^{2}\geq F(\hat{\beta}_{(X,y)}(\lambda))\geq\lambda_{1}||{\hat{\beta}_{(X,y)}(\lambda)}||_{1}+\lambda_{2}||{\hat{\beta}_{(X,y)}(\lambda)}||_{2}^{2}.$$

Therefore, for any problem instance $P$, the solution of the training optimization problem $\hat{\beta}_{(X,y)}(\lambda)$ has bounded norm, i.e. $||{\hat{\beta}_{(X,y)}(\lambda)}||_{1},||{\hat{\beta}_{(X,y)}(\lambda)}||_{2}^{2}\leq\frac{R^{2}}{2\lambda_{\min}}$, which implies

$$h(\lambda,P)=\frac{1}{2m}||{y_{\text{val}}-\hat{\beta}_{(X,y)}(\lambda)X_{\text{val}}}||_{2}^{2}\leq\frac{1}{2m}||{y_{\text{val}}}||_{2}^{2}+\frac{1}{2m}||{\hat{\beta}_{(X,y)}(\lambda)}{X_{\text{val}}}||_{2}^{2}\leq H,$$

for some constant $H$ (depends only on $R$ and $\lambda_{\min}$). ∎

C.2 Tuning the ElasticNet – Online learning

At a high level, the plan is to show dispersion (Definition 4) using the general recipe developed in [BDP20]. The recipe may be summarized at a high level as follows.

• Bound the probability density of the random set of discontinuities of the loss functions. Intuitively this corresponds to computing the average number of loss functions that may be discontinuous along a path connecting any two points within distance $\epsilon$ in the domain.

• Use a VC-dimension based uniform convergence argument to transform this into a bound on the dispersion of the loss functions.

Formally, we have the following theorems from [BDP20], which show how to use this technique when the discontinuities are roots of a random polynomial with bounded coefficients. The theorems implement steps S1 and S2 of the above recipe respectively.

Theorem C.7 ([BDP20]).

Consider a random degree $d$ polynomial $\phi$ with leading coefficient 1 and subsequent coefficients which are real of absolute value at most $R$, whose joint density is at most $\kappa$. There is an absolute constant $K_{0}$ depending only on $d$ and $R$ such that every interval $I$ of length $\leq\epsilon$ satisfies Pr($\phi$ has a root in $I$) $\leq\kappa\epsilon/K_{0}$.

Theorem C.8 ([BDP20]).

Let $l_{1},\dots,l_{T}:\mathbb{R}\rightarrow\mathbb{R}$ be independent piecewise $L$-Lipschitz functions, each having at most $K$ discontinuities. Let $D(T,\epsilon,\rho)=|\{1\leq t\leq T\mid l_{t}\text{ is not }L\text{-Lipschitz on }[\rho-\epsilon,\rho+\epsilon]\}|$ be the number of functions that are not $L$-Lipschitz on the ball $[\rho-\epsilon,\rho+\epsilon]$. Then we have $E[\max_{\rho\in\mathbb{R}}D(T,\epsilon,\rho)]\leq\max_{\rho\in\mathbb{R}}E[D(T,\epsilon,\rho)]+O(\sqrt{T\log(TK)})$.

The following lemma provides useful extension to Lemma 2.1 for our online learning results.

Lemma C.9.

Let $A$ be an $r\times s$ matrix with $R$-bounded max-norm, i.e. $||A||_{\infty,\infty}=\max_{i,j}|A_{ij}|\leq R$. Then each entry of the matrix $(A^{T}A+\lambda I_{s})^{-1}$ is a rational polynomial $P_{ij}(\lambda)/Q(\lambda)$ for $i,j\in[s]$ with each $P_{ij}$ of degree at most $s-1$, $Q$ of degree $s$, and all the coefficients have absolute value at most $r^{s}(Rs)^{2s}$.

Proof.

Let $G=A^{T}A$ be the Gram matrix.

$$|G_{ij}|=|\sum_{k}{A_{ki}A_{kj}}|\leq\sum_{k}{|A_{ki}A_{kj}|}\leq rR^{2},$$

by the triangle inequality and using $\max_{i,j}|A_{ij}|\leq R$. The determinant $\textsc{det}(A^{T}A+\lambda I_{s})$ is a sum of $s!\leq s^{s}$ signed terms, each a product of $s$ elements of the form $G_{ij}$ or $G_{ii}+\lambda$. Thus, in each of the $s!$ terms, the coefficient of $\lambda^{k}$ is a sum of at most ${s\choose s-k}\leq s^{k}\leq s^{s}$ expressions of the form $\Pi_{(i,j)\in S}G_{ij}$ with $|S|\leq s-k$. Now $|\Pi_{(i,j)\in S}G_{ij}|\leq(rR^{2})^{|S|}\leq(rR^{2})^{s}$, and by triangle inequality the coefficient of $\lambda^{k}$ is upper bounded by $(rR^{2})^{s}\cdot s^{k}\cdot s^{s}$ for any $k$. This establishes the bound on the coefficients of $Q(\lambda)$. A similar argument implies the upper bound for each $P_{ij}(\lambda)$. ∎

We will also need the following result, which is a simple extension of Lemma 24 from [BS21].

Lemma C.10.

Suppose $X$ and $Y$ are real-valued random variables taking values in $[m,m+M]$ for some $m,M\in\mathbb{R}^{+}$ and suppose that their joint distribution is $\kappa$-bounded. Let $c$ be an absolute constant. Then,

• $Z=X+Y$ is drawn from a $K_{1}\kappa$-bounded distribution, where $K_{1}\leq M$.

• $Z=XY$ is drawn from a $K_{2}\kappa$-bounded distribution, where $K_{2}\leq M/m$.

• $Z=X-Y$ is drawn from a $K_{1}\kappa$-bounded distribution, where $K_{1}\leq M$.

• $Z=X+c$ has a $\kappa$-bounded distribution, and $Z=cX$ has a $\frac{\kappa}{|c|}$-bounded distribution.

Proof.

Let $f_{X,Y}(x,y)$ denote the joint density of $X,Y$. (i) and (ii) are immediate from Lemma 24 from [BS21], (iii) is a simple extension. Indeed, the cumulative density function for $Z$ is given by

$$F_{Z}(z) =\text{Pr}(Z\leq z)=\text{Pr}(X-Y\leq z)=\text{Pr}(X\leq z+Y)$$

$$=\int_{m}^{m+M}\int_{m}^{z+y}f_{X,Y}(x,y)dxdy.$$

The density function for $Z$ can be obtained using Leibniz’s rule as

$$f_{Z}(z)=\frac{d}{dz}F_{Z}(z) =\frac{d}{dz}\int_{m}^{m+M}\int_{m}^{z+y}f_{X,Y}(x,y)dxdy$$

$$=\int_{m}^{m+M}\left(\frac{d}{dz}\int_{m}^{y}f_{X,Y}(x,y)dx+\frac{d}{dz}\int_{0}^{z}f_{X,Y}(t+y,y)dt\right)dy$$

$$=\int_{m}^{m+M}f_{X,Y}(z+y,y)dy$$

$$\leq\int_{m}^{m+M}\kappa dy$$

$$=M\kappa.$$

Finally, (iv) follows from simple change of variable manipulations (e.g. Theorem 22 of [BDP20]). ∎

Theorem 3.3 (restated). Assume that the predicted variable and all feature values are bounded by an absolute constant $R$, i.e. $\max\{||X^{(i)}||_{\infty,\infty},||y^{(i)}||_{\infty},||X_{\text{val}}^{(i)}||_{\infty,\infty},||y_{\text{val}}^{(i)}||_{\infty}\}\leq R$. Suppose the predicted variables $y^{(i)}$ in the training set are drawn from a joint $\kappa$-bounded distribution. Let $l_{1},\dots,l_{T}:(0,\lambda_{\max})^{2}\rightarrow\mathbb{R}_{\geq 0}$ denote an independent sequence of losses (e.g. fresh randomness is used to generate the validation set features in each round) as a function of the ElasticNet regularization parameter $\lambda=(\lambda_{1},\lambda_{2})$, $l_{i}(\lambda)=l_{r}(\hat{\beta}^{(X^{(i)},y^{(i)})}_{\lambda,f_{EN}},(X_{\text{val}}^{(i)},y_{\text{val}}^{(i)}))$. The sequence of functions is $\frac{1}{2}$-dispersed, and there is an online algorithm with $\tilde{O}(\sqrt{T})$ expected regret. The result also holds for loss functions adjusted by information criteria AIC and BIC.

Proof.

We start with the piecewise-decomposable characterization of the dual class function in Theorem 2.2. On any fixed problem instance $P\in\Pi_{m,n}$, as the parameter $\lambda$ is varied in the loss function $\ell_{\text{EN}}(\cdot,P)$ of ElasticNet trained with regularization parameter $\lambda=(\lambda_{1},\lambda_{2})$, we have the following piecewise structure. There are $k=p3^{p}$ boundary functions $g_{1},\dots,g_{k}$ for which the transition boundaries are algebraic curves $r_{i}(\lambda_{1},\lambda_{2})$, where $r_{i}$ is a polynomial with degree 1 in $\lambda_{1}$ and at most $p$ in $\lambda_{2}$. Also the piece function $f_{\mathbf{b}}$ for each sign pattern $\mathbf{b}\in\{0,1\}^{k}$ is a rational polynomial function $\frac{q_{1}^{b}(\lambda_{1},\lambda_{2})}{q_{2}^{b}(\lambda_{2})}$, where $q_{1}^{\mathbf{b}},q_{2}^{\mathbf{b}}$ have degrees at most $2p$, and corresponds to a fixed signed equicorrelation set $(\mathcal{E},s)$. To show online learnability, we will examine this piecewise structure more closely – in particular analyse how the structure varies when the predicted variable is drawn from a smooth distribution.

In order to show dispersion for the loss functions $\{l_{i}(\lambda)\}$, we will use the recipe of [BDP20] and bound the worst rate of discontinuities between any pair of points $\lambda=(\lambda_{1},\lambda_{2})$ and $\lambda^{\prime}=(\lambda^{\prime}_{1},\lambda^{\prime}_{2})$ with $||\lambda-\lambda^{\prime}||_{2}\leq\epsilon$ along the axis-aligned path $\lambda\rightarrow(\lambda^{\prime}_{1},\lambda_{2})\rightarrow\lambda^{\prime}$. First observe that the only possible points at which $l_{i}(\lambda)$ may be discontinuous are

• $(\lambda_{1},\lambda_{2})$ such that $r_{i}(\lambda_{1},\lambda_{2})=0$ corresponding to some boundary function $g_{i}$.

• $(\lambda_{1},\lambda_{2})$ such that $q_{2}^{\mathbf{b}}(\lambda_{2})=0$ corresponding to some piece function $f_{\mathbf{b}}$.

Fortunately the discontinuity of type (b) does not occur for $\lambda_{2}>0$. From the ElasticNet characterization in Lemma C.1, and using Lemma 2.1, we know that $q_{2}(\lambda_{2})=\Pi_{j\in[|\mathcal{E}|]}(\Lambda_{j}+\lambda_{2})$, where $(\Lambda_{j})_{j\in[|\mathcal{E}|]}$ are non-negative eigenvalues of the positive semi-definite matrix ${X^{(i)}_{\mathcal{E}}}^{T}X^{(i)}_{\mathcal{E}}$. It follows that $q_{2}^{\mathbf{b}}$ does not have positive zeros (for any sign vector $\mathbf{b}$).

Therefore it suffices to locate boundaries of type (a). To this end, we have two subtypes corresponding to a variable entering or leaving the equicorrelation set.

Addition of $j\notin\mathcal{E}$. As observed in the proof of Theorem 2.2, a variate $j\notin\mathcal{E}$ can enter the equicorrelation set $\mathcal{E}$ only for $(\lambda_{1},\lambda_{2})$ satisfying

$$\lambda_{1}=\frac{\bm{x}_{j}^{T}X_{\mathcal{E}}({X_{\mathcal{E}}}^{T}X_{\mathcal{E}}+\lambda_{2}I_{|\mathcal{E}|})^{-1}{X_{\mathcal{E}}}^{T}y-\bm{x}_{j}^{T}y}{\bm{x}_{j}^{T}X_{\mathcal{E}}({X_{\mathcal{E}}}^{T}X_{\mathcal{E}}+\lambda_{2}I_{|\mathcal{E}|})^{-1}s\pm 1}.$$

For fixed $\lambda_{2}$, the distribution of $\lambda_{1}$ at which the discontinuity occurs for insertion of $j$ is $K_{1}\kappa$-bounded (by Lemma C.10) for some constant $K_{1}$ that only depends on $R,m,p$ and $\lambda_{\max}$. This implies an upper bound of $K_{1}\kappa\epsilon$ on the expected number of discontinuities corresponding to $j$ along the segment $\lambda\rightarrow(\lambda^{\prime}_{1},\lambda_{2})$ for any $j,\mathcal{E}$.

For constant $\lambda_{1}$, we can use Lemma 2.1 and a standard change of variable argument (e.g. Theorem 22 of [BDP20]) to conclude that the discontinuties lie at the roots of a random polynomial in $\lambda_{2}$ of degree $|\mathcal{E}|$, leading coefficient $1$, and bounded random coefficients with $K_{2}\kappa$-bounded density for some constant $K_{2}$ (that only depends on $R,m,p$ and $\lambda_{\max}$). By Theorem C.7, the expected number of discontinuities along the segment $(\lambda^{\prime}_{1},\lambda_{2})\rightarrow\lambda^{\prime}$ is upper bounded by $K_{2}K_{p}\kappa\epsilon$ ($K_{p}$ only depends on $p$). This implies that the expected number of Lipschitz violations between $\lambda$ and $\lambda^{\prime}$ along the axis aligned path is $\tilde{O}(\kappa\epsilon)$ and completes the first step of the recipe in this case ($\tilde{O}$ notation suppresses terms in $R,m,p$ and $\lambda_{\max}$ as constants).

Removal of $j^{\prime}\in\mathcal{E}$. The second case, when a variate $j^{\prime}\in\mathcal{E}$ leaves the equicorrelation set $\mathcal{E}$ for $(\lambda_{1},\lambda_{2})$ satisfying

$$\lambda_{1}(({X_{\mathcal{E}}}^{T}X_{\mathcal{E}}+\lambda_{2}I_{|\mathcal{E}|})^{-1}s)_{j^{\prime}}=(({X_{\mathcal{E}}}^{T}X_{\mathcal{E}}+\lambda_{2}I_{|\mathcal{E}|})^{-1}{X_{\mathcal{E}}}^{T}y)_{j^{\prime}},$$

also yields the same bound using the above arguments. Putting together, and noting that we have at most $p3^{p}$ distinct curves each with $\tilde{O}(\kappa\epsilon)$ expected number of intersections with the axis aligned path $\lambda\rightarrow\lambda^{\prime}$, the total expected number of discontinuities is also $\tilde{O}(\kappa\epsilon)$. This completes the first step (S1) of the above recipe.

We use Theorem 9 of [BDP20] to complete the second step of the recipe, which employs a VC-dimension argument for $K^{\prime}$ algebraic curves of bounded degrees (here degree is at most $p+1$) to conclude that the expected worst number of discontinuties along any axis-aligned path between any pair of points $\leq\epsilon$ apart is at most $\tilde{O}(\epsilon T)+O(\sqrt{T\log K^{\prime}T})$. $K^{\prime}\leq p3^{p}$ as shown above. This implies that the sequence of loss functions is $\frac{1}{2}$-dispersed, and further there is an algorithm (Algorithm 4 of [BDV18]) that achieves $\tilde{O}(\sqrt{T})$ expected regret.

Finally note that loss functions with AIC and BIC have the same dual class piecewise structure, and therefore the above analysis applies. The only difference is that the value of the piece functions $f_{\mathbf{b}}$ are changed by a constant (in $\lambda$), $K_{m,p}\leq p\log m$. The piece boundaries are the same, and are therefore $\frac{1}{2}$-dispersed as above. The range of the loss functions is now $[0,K_{m,p}+1]$, so the same algorithm (Algorithm 4 of [BDV18]) again achieves $\tilde{O}(\sqrt{T})$ expected regret. ∎

Appendix D Lemmas and proof details for Section 4

We will first extend the structure for the ElasticNet regression loss functions shown in Theorem 2.2 to the classification setting. The main new challenge is that there are additional discontinuties due to thresholding the loss function needed for binary classification, which intuitively makes the loss more jumpy and discontinuous as a function of the regularization parameters.

Lemma D.1.

Let $\mathcal{L}$ be a set of functions $\{l_{\lambda,\tau}:\Pi_{m,p}\rightarrow\mathbb{R}_{\geq 0}\mid\lambda\in\mathbb{R}^{+}\times\mathbb{R}_{\geq 0},\tau\in\mathbb{R}\}$ that map a regression problem instance $P\in\Pi_{m,p}$ to the validation classification loss $\ell_{\text{EN}}^{c}(\lambda,P,\tau)$ of ElasticNet trained with regularization parameter $\lambda=(\lambda_{1},\lambda_{2})$ and threshold parameter $\tau$. The dual class $\mathcal{L}^{*}$ is $(\mathcal{F},\mathcal{G},(m+p)3^{p})$-piecewise decomposable, with $\mathcal{F}=\{f_{c}:\mathcal{L}\rightarrow\mathbb{R}\}$ consisting of constant functions $f_{c}:l_{\lambda,\tau}\mapsto c$, where $c\in\mathbb{R}_{\geq 0}$, and $\mathcal{G}=\{g_{r}:\mathcal{L}\rightarrow\{0,1\}\}$ consisting of semi-algebraic sets bounded by algebraic varieties $g_{r}:l_{\lambda,\tau}\mapsto\mathbb{I}\{r(\lambda_{1},\lambda_{2},\tau)<0\}$, where $r$ is a polynomial of degree 1 in $\lambda_{1}$ and $\tau$, and at most $p$ in $\lambda_{2}$.

Proof.

By Lemma C.1, the EN coefficients $\hat{\beta}_{EN}$ are fixed given the signed equicorrelation set $\mathcal{E},s$. As in Theorem 2.2, we have $\leq p3^{p}$ boundaries $\mathcal{G}_{1}$ corresponding to a change in the equicorrelation set, but the value of the loss also changes when a prediction vector coefficient $\mu_{j}=(X_{\text{val}})_{j}\hat{\beta}_{EN}$ cross the threshold $\tau$. This is given by $(c_{1}-c_{2}{\lambda_{1}})_{j}=\tau$ where

$$c_{1}=(X_{\text{val}})_{\mathcal{E}}(X_{\mathcal{E}}^{T}X_{\mathcal{E}}+\lambda_{2}I_{|\mathcal{E}|})^{-1}X_{\mathcal{E}}^{T}y,$$

and

$$c_{2}=(X_{\text{val}})_{\mathcal{E}}(X_{\mathcal{E}}^{T}X_{\mathcal{E}}+\lambda_{2}I_{|\mathcal{E}|})^{-1}s.$$

Therefore, $\hat{\mu}_{j}=0$ corresponds to the 0-set of $(c_{1}-c_{2}{\lambda_{1}})_{j}-\tau$. By an application of Lemma 2.1, this is an algebraic variety with degree at most $|\mathcal{E}|$ in $\lambda_{2}$ and degree $1$ in $\lambda_{1}$ and $\tau$. There are at most $m3^{p}$ such boundary functions $\mathcal{G}_{2}$, corresponding to all possibilities of $j,\mathcal{E},s$. For a point $(\lambda_{1},\lambda_{2},\tau)$ with a fixed sign pattern of boundary functions in $\mathcal{G}_{1}\cup\mathcal{G}_{2}$, the EN coefficients are fixed and also all the predictions on the validation set are fixed. The classification loss

$$l_{c}(\hat{\beta}_{\lambda,f},(X,y),\tau)=\frac{1}{m}\sum_{i=1}^{m}(y_{i}-\textup{{sgn}}(\langle X_{i},\hat{\beta}_{\lambda,f}\rangle-\tau))^{2}$$

is therefore constant in each piece. Applying Theorem 2.2 shows the claimed stucture for the ElasticNet classification (dual class) loss function. ∎

The above piecewise decomposable structure is helpful in bounding the pseudodimension for the ElasticNet based classifier. For the special cases of Ridge and LASSO we obtain the pseudodimension bounds from first principles.

Theorem 4.1 (restated). Let $\mathcal{H}_{\text{Ridge}}^{c}$, $\mathcal{H}_{\text{LASSO}}^{c}$ and $\mathcal{H}_{\text{EN}}^{c}$ denote the set of loss functions for classification problems with at most $m$ examples and $p$ features, with Ridge, LASSO and ElasticNet regularization respectively.

• $\textsc{Pdim}(\mathcal{H}_{\text{Ridge}}^{c})=O(\log mp)$

• $\textsc{Pdim}(\mathcal{H}_{\text{LASSO}}^{c})=O(p\log m)$. Further, in the overparameterized regime ($p\gg m$), we have that $\textsc{Pdim}(\mathcal{H}_{\text{LASSO}}^{c})=O(m\log\frac{p}{m})$.

• $\textsc{Pdim}(\mathcal{H}_{\text{EN}}^{c})=O(p^{2}+p\log m)$.

Proof.

For Ridge regression, the estimator $\hat{\beta}_{\lambda,f_{2}}$ on the dataset $(X^{(i)},y^{(i)})$ is given by the following closed form

$$\hat{\beta}_{\lambda,f_{2}}=({X^{(i)}}^{T}X^{(i)}+\lambda I_{p_{i}})^{-1}{X^{(i)}}^{T}y^{(i)},$$

where $I_{p_{i}}$ is the ${p_{i}}\times{p_{i}}$ identity matrix. By Lemma 2.1 each coefficient $(\hat{\beta}_{\lambda,f_{2}})_{k}$ of the estimator $\hat{\beta}_{\lambda,f_{2}}$ is a rational polynomial in $\lambda$ of the form $P_{k}(\lambda)/Q(\lambda)$, where $P_{k},Q$ are polynomials of degrees at most $p_{i}-1$ and $p_{i}$ respectively. Thus the prediction on any example $(X_{\text{val}}^{(i)})_{j}$ in the validation set $(X_{\text{val}}^{(i)},y_{\text{val}}^{(i)})$ of any problem instance $P^{(i)}$ can change at most $p_{i}\leq p$ times as $\lambda$ is varied. Recall there are $m_{i}^{\prime}\leq m$ examples in any validation set. This implies we have at most $mnp$ distinct values of the loss function over the $n$ problem instances. The pseudo-dimension $n$ therefore satisfies $2^{n}\leq mnp$, or $n=O(\log mp)$.

Prior work [EHJT04] shows that the optimal vector $\hat{\beta}\in\mathbb{R}^{p}$ evolves piecewise linearly with $\lambda$, i.e. $\exists\lambda^{(0)}=0<\lambda^{(1)}<\dots<\lambda^{(q)}=\infty$ and $\gamma_{0},\gamma_{1},\dots,\gamma_{q-1}\in\mathbb{R}^{p}$ such that

$$\hat{\beta}_{\lambda,f_{1}}=\hat{\beta}_{\lambda^{(k)},f_{1}}+(\lambda-\lambda^{(k)})\gamma_{k}$$

for $\lambda^{(k)}\leq\lambda\leq\lambda^{(k+1)}$. Each piece corresponds to the addition or removal of at least one of $p$ coordinates to the active set of covariates with maximum correlation. For any data point $x_{j}$, $1\leq j\leq m$, and any piece $[\lambda^{(k)},\lambda^{(k+1)})$ we have that $x_{j}\hat{\beta}$ is monotonic since $\hat{\beta}$ varies along a fixed vector $\gamma_{k}$, and therefore can have at most one value of $\lambda$ where the predicted label $\hat{y}$ changes. This gives an upper bound of $mq$ on the total number of discontinuities on any single problem instance $(X^{(i)},y^{(i)},X_{\text{val}}^{(i)},y_{\text{val}}^{(i)})$, where $q$ is the number of pieces in the solution path. By Lemma 6 of [Tib13], we have the number pieces in the solution path $q\leq 3^{p}$. Also for the overparameterized regime $p\gg m$, we have the property that there are at most $m-1$ variables in the active set for the entire sequence of solution paths (Section 7, [EHJT04]). Thus, we have that $q\leq m{p\choose m-1}\leq(\frac{ep}{m})^{m}$ in this case.

Over $n$ problem instances, the pseudo-dimension satisfies $2^{n}\leq mqn$, or $n=O(\log mq)$. Substituting the above inequalities for $q$ completes the proof.

The proof of Theorem follows the same arguments as Theorem 3.1, using Lemma D.1 instead of Theorem 2.2. By Lemma D.1, the dual class $\mathcal{H}_{\text{EN}}^{*}$ of $\mathcal{H}_{\text{EN}}$ is $(\mathcal{F},\mathcal{G},(p+m)3^{p})$-piecewise decomposable, with $\mathcal{F}=\{f_{c}:\mathcal{L}\rightarrow\mathbb{R}\}$ consisting of constant functions $f_{c}:l_{\lambda}\mapsto c$, where $c\in\mathbb{R}_{\geq 0}$, and $\mathcal{G}=\{g_{r}:\mathcal{L}\rightarrow\{0,1\}\}$ consisting of polynomial thresholds $g_{r}:u_{\lambda}\mapsto\mathbb{I}\{r(\lambda_{1},\lambda_{2},\tau)<0\}$, where $r$ is a polynomial of degree 1 in $\lambda_{1}$ and $\tau$, and at most $p$ in $\lambda_{2}$.

Now it is easy to see that $\textsc{Pdim}(\mathcal{F}^{*})=O(1)$ (in particular, a consequence of Lemma C.4), and by Lemma C.5 the VC dimension of the dual boundary class is $d_{\mathcal{G}^{*}}=O(p)$. A straightforward application of Theorem C.2 yields

$$\textsc{Pdim}(\mathcal{H})=O(p\log p+p\log((p+m)3^{p}))=O(p^{2}+p\log m).$$

∎

We will now restate and prove Theorem 4.2. This implies that under smoothness assumptions on the data distribution we can learn the data-dependent optimal regularization parameter in the online setting.

Theorem 4.2 (restated). Suppose Assumptions 1 and 3 hold. Let $l_{1},\dots,l_{T}:(0,H]^{d}\times[-H,H]\rightarrow\mathbb{R}$ denote an independent sequence of losses as a function of the regularization parameter $\lambda$, $l_{i}(\lambda,\tau)=l_{c}(\hat{\beta}_{\lambda,f},(X^{(i)},y^{(i)}),\tau)$. The sequence of functions is $\frac{1}{2}$-dispersed, and there is an online algorithm with $\tilde{O}(\sqrt{T})$ expected regret, if $f$ is given by

• $f=f_{1}$ (LASSO),

• $f=f_{2}$ (Ridge), or

• $f=f_{EN}$ (ElasticNet).

Proof.

On any dataset $(X,y)$, the predictions are given by the coefficients of the prediction vector $\hat{\mu}=X(X^{T}X+\lambda I_{p})^{-1}X^{T}y$. Note that by Lemma 2.1 $(X^{T}X+\lambda I_{p})^{-1}$, and therefore $X_{\text{val}}(X^{T}X+\lambda I_{p})^{-1}X^{T}y$, has each element of the form $P_{j}(\lambda)/Q(\lambda)$ with degree of each $P_{j}$ at most $p-1$ and degree of $Q$ at most $p$. Further, for a fixed $\tau$, by using Lemma C.9 and a change of variables, we have that $\hat{\mu}_{j}=\tau$ is polynomial equation in $\lambda$ with degree $p$ with bounded coefficients that have $K\kappa$-bounded density for some constant $K$ (that depends on $R$, $H$, $m$ and $p$, but not on $\kappa$) and leading coefficient 1. Further, for fixed $\lambda\in(0,H]$, by Lemma C.10, the discontinuities over $\tau$ again have $\tilde{O}(\kappa)$-bounded density. This completes step S1 of the recipe from [BDP20], and Theorem C.7 gives a bound on the expected number of discontinuities in any interval $I$ (over $\lambda$).

To complete step S2, note that the loss function $l_{i}$ on any instance has at most $pm$ discontinuities, since each coefficient $\hat{\mu}_{j}$ of the prediction vector can change sign at most $p$ times as $\lambda$ is varied. This implies the VC-dimension argument (Theorem C.8) applies and the expected maximum number of discontinuities in any interval of width $\epsilon$ is $O(\epsilon T)+O(\sqrt{T\log(mpT)})$, which is $\tilde{O}(\epsilon T)$ for $\epsilon\geq 1/\sqrt{T}$. Thus, using the recipe from [BDP20], we have shown that the sequence of loss functions is $\frac{1}{2}$-dispersed. This further implies that Algorithm 1, which implements the Continuous Exp-Weights algorithm of [BDV18] for setting the regularization parameter, achieves $\tilde{O}(\sqrt{T})$ expected regret ([BDV18], Theorem 1).

Since the data-distribution is in particular assumed to be continuous, by Lemma 4 of [Tib13] we know that the LASSO solutions are unique with probability 1. Moreover if $\mathcal{E}\subseteq[p]$ denotes the equicorrelation set of variables (i.e. covariates with the maximum absolute value of correlation), and $s\in\{-1,1\}^{|\mathcal{E}|}$ the sign vector (i.e. the sign of the correlations of the covariates in $\mathcal{E}$), then the LASSO prediction vector $\hat{\mu}=X_{\text{val}}\hat{\beta}$ is a linear function of regularization parameter $\lambda$ given by

$$\hat{\mu}=c_{1}-c_{2}\lambda,$$

where $c_{1}=(X_{\text{val}})_{\mathcal{E}}(X_{\mathcal{E}}^{T}X_{\mathcal{E}})^{-1}X_{\mathcal{E}}^{T}y$ and $c_{2}=(X_{\text{val}})_{\mathcal{E}}(X_{\mathcal{E}}^{T}X_{\mathcal{E}})^{-1}s$. Thus for any fixed $\mathcal{E},s$ (corresponding to a unique piece in the solution path for LARS-LASSO), we have at most one discontinuity corresponding to $\hat{\mu}_{j}=\tau$, and the location of this discontinuity has a $K\kappa$-bounded distribution (for constant $K$ independent of $\kappa$) by an application of Lemma C.10. Thus, the probability that this discontinuity is located along some axis aligned path $I$ of length $\epsilon$ is at most $K\kappa\epsilon$. A union bound over $j\in[m]$, and over $3^{p}$ choices of $\mathcal{E},s$ (for example, Lemma 6 in [Tib13]) gives the probability of a discontinuity along $I$ is at most $m3^{p}K\kappa\epsilon$. This completes step S1 of the recipe above.

Now each loss function $l_{i}$ has at most $m3^{p}$ discontinuities, and therefore by a VC-dimension argument (Theorem 9 of [BDP20]), the expected maximum number of discontinuities along any axis-aligned path of total length $\epsilon$ is $\tilde{O}(\epsilon T)+O(\sqrt{T(p+\log(mT))})$, which is $\tilde{O}(\epsilon T)$ for $\epsilon\geq 1/\sqrt{T}$. This completes step S2 of the recipe from [BDP20], and we have shown that the sequence of loss functions is $\frac{1}{2}$-dispersed. As in Theorem 4.2, this implies that Algorithm 1 achieves $\tilde{O}(\sqrt{T})$ expected regret ([BDV18], Theorem 1).

While we use the worst case bound on the number of solution paths here, algorithmically we can use LARS-LASSO on the given dataset, which is much faster in practice and the running time scales linearly with the actual number of solution paths $q$ (typically $q\ll 3^{p}$).

The proof uses the piecewise decomposable structure proved in Lemma D.1, and establishes dispersion using joint smoothness of $X_{\text{val}}^{(i)}$ instead of $y^{(i)}$ (as in the proof of Theorem 4.2 (i)). The recipe from [BDP20] can be used along a 3D axis-aligned path from $(\lambda_{1},\lambda_{2},\tau)\rightarrow(\lambda_{1}^{\prime},\lambda_{2}^{\prime},\tau^{\prime})$. Lemma C.10 may be used to show bounded density of discontinuities along $\tau$ (keeping $\lambda_{1},\lambda_{2}$ fixed). To complete step S2 of the recipe we can use Theorem 7 from [BS21]. The arguments are otherwise very similar to the proof of Theorem 3.3, and are omitted for brevity. ∎