# Multivariate Convex Regression at Scale

Wenyu Chen MIT Operations Research Center. email<EMAIL_ADDRESS>Massachusetts
Institute of Technology  Rahul Mazumder MIT Sloan School of Management,
Operations Research Center and MIT Center for Statistics. email:
<EMAIL_ADDRESS>Massachusetts Institute of Technology

(May, 2020)

###### Abstract

We present new large-scale algorithms for fitting a subgradient regularized
multivariate convex regression function to $n$ samples in $d$ dimensions—a key
problem in shape constrained nonparametric regression with widespread
applications in statistics, engineering and the applied sciences. The
infinite-dimensional learning task can be expressed via a convex quadratic
program (QP) with $O(nd)$ decision variables and $O(n^{2})$ constraints. While
instances with $n$ in the lower thousands can be addressed with current
algorithms within reasonable runtimes, solving larger problems (e.g.,
$n\approx 10^{4}$ or $10^{5}$) is computationally challenging. To this end, we
present an active set type algorithm on the dual QP. For computational
scalability, we perform approximate optimization of the reduced sub-problems;
and propose randomized augmentation rules for expanding the active set.
Although the dual is not strongly convex, we present a novel linear
convergence rate of our algorithm on the dual. We demonstrate that our
framework can approximately solve instances of the convex regression problem
with $n=10^{5}$ and $d=10$ within minutes; and offers significant
computational gains compared to earlier approaches.

##  1 Introduction

We consider the task of computing a multivariate convex regression function
given samples $(y_{i},\bm{x}_{i})$, $i\in[n]:=\\{1,\ldots,n\\}$ with response
$y_{i}\in\mathbb{R}$ and covariates $\bm{x}_{i}\in\mathbb{R}^{d}$. The convex
least squares fit [35, 21] is given by the following infinite-dimensional
problem

$\hat{\varphi}=\arg\min_{\varphi\in\mathcal{F}}~{}\frac{1}{n}\sum_{i=1}^{n}(y_{i}-\varphi(\bm{x}_{i}))^{2}$
(1)

where, $\mathcal{F}$ is the set of all real-valued convex functions in
$\mathbb{R}^{d}$, $\varphi(\bm{x}_{i})$ is the value of the function $\varphi$
at point $\bm{x}_{i}$. Criterion (1) seeks to find a convex function
$\varphi(\bm{x})$, with the best data-fidelity in terms of the squared error
loss on the observed samples $(y_{i},\bm{x}_{i}),i\in[n]$. Different from
common nonparametric function estimation methods under smoothness constraints
[38, 13], Problem (1) is an instance of shape constrained nonparametric
regression [33]—a subfield of nonparametric statistics where learning is
carried out under qualitative shape constraints on the underlying function.
This topic has received significant attention over the recent years—see for
example, the special issue in Statistical Science [34] for a nice overview.

Convex regression is widely used in economics, operations research,
statistical learning and engineering applications. In economics applications,
for example, convexity/concavity arise in modeling utility and production
functions, consumer preferences [37, 18], among others. In some stochastic
optimization problems, value functions are taken to be convex [36]. See also
the works of [16, 39, 3] for other important applications of convex
regression.

Statistical methodology underlying convex regression has been investigated by
several authors [35, 21, 3, 14, 15, 20]. More recently, to realize the full
potential of such methods in statistics/data science practice there is a
growing interest in developing efficient algorithms for this problem–see for
e.g. [25, 1, 6, 22]—the focus of the current paper is to further advance the
computational frontiers of convex regression.

An appealing aspect of the functional optimization Problem (1) is that it can
be reduced [35, 25] to a finite dimensional convex quadratic program (QP):

$\displaystyle\underset{\phi_{1},\ldots,\phi_{n};\bm{\xi}_{1},\ldots,\bm{\xi}_{n}}{\text{minimize}}$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}(y_{i}-\phi_{i})^{2}$ (2) s.t.
$\displaystyle\phi_{j}-\phi_{i}\geq\langle\bm{x}_{j}-\bm{x}_{i},\bm{\xi}_{i}\rangle,\quad\forall(i,j)\in\Omega$

where $\phi_{1},\ldots,\phi_{n}\in\mathbb{R}$,
$\bm{\xi}_{1},\ldots,\bm{\xi}_{n}\in\mathbb{R}^{d}$, and
$\Omega:=\\{(i,j):1\leq i,j\leq n,i\neq j\\}$. In (2),
$\phi_{i}=\varphi(\bm{x}_{i})$ and $\bm{\xi}_{i}$ is a subgradient of
$\varphi(\bm{x})$ at $\bm{x}=\bm{x}_{i}$. Problem (2) has $O(nd)$ variables
and $O(n^{2})$ constraints and becomes computationally challenging when $n$ is
large. For the convex regression problem to be statistically meaningful [20,
25], we consider cases with $n\gg d$ (with $d\sim\log n$). Off-the-shelf
interior point methods [35] for (2) are limited to instances where $n$ is at
most a few hundred. [1] consider a regularized version of (2) (i.e.,
($P_{0}$), below) and propose parallel algorithms to solve instances with
$n\approx 1,600$ leveraging commercial QP solvers (like Mosek). [25] use an
alternating direction method of multipliers (ADMM) [8]-based algorithm that
can address problems up to $n\approx 3,000$. Recently, [22] propose a
different ADMM method and also a proximal augmented Lagrangian method where
the subproblems are solved by the semismooth Newton method—they address
instances with $n\approx 3,000$. Algorithms based on nonconvex optimization
have been proposed to learn convex functions that are representable as a
piecewise maximum of $k$-many hyperplanes [16, 3, 12]—these are interesting
approaches, but they may not lead to an optimal solution for the convex
regression convex optimization problem. Recently, [6] present a cutting plane
or constraint generation-type algorithm for (2): At every iteration, they
solve a reduced QP by considering a subset of constraints. Leveraging
capabilities of commercial solvers (e.g., Gurobi), this can approximately
solve instances of (2) with $n\approx 10^{4}$-$10^{5}$. In this paper, we also
present an active-set type method, but our approach differs from [6], as we
discuss below. We also establish novel computational guarantees for our
proposed approach.

The convex least squares estimator (2) may suffer in terms of out-of-sample
performance due to spurious boundary effects [25, 2]. This problem can be
alleviated in part with additional subgradient regularization. To this end, we
consider a subgradient regularized-version of (2) given by:

$\displaystyle\begin{aligned}
&\underset{\phi_{1},\ldots,\phi_{n};\bm{\xi}_{1},\ldots,\bm{\xi}_{n}}{\text{minimize}}&&\frac{1}{n}\sum_{i=1}^{n}(y_{i}-\phi_{i})^{2}+\frac{\rho}{n}\sum_{i=1}^{n}\lVert\bm{\xi}_{i}\rVert^{2}\\\
&\quad\,\text{s.t.}&&\phi_{j}-\phi_{i}\geq\langle\bm{x}_{j}-\bm{x}_{i},\bm{\xi}_{i}\rangle,\quad\forall(i,j)\in\Omega\end{aligned}$
($P_{0}$)

where, we impose an additional $\ell_{2}$-based regularization on the
subgradients $\\{\bm{\xi}_{i}\\}_{1}^{n}$ of $\varphi$; with $\rho>0$ being
the regularization parameter. Formulation ($P_{0}$) also appears in [1].
Estimation error properties of a subgradient regularized convex least squares
estimator have been shown in [25]. As seen in our numerical experiments, a
choice of $\rho>0$ can lead to superior out-of-sample performance compared to
$\rho=0$.

Our approach: In this paper, we focus on solving ($P_{0}$). Problem ($P_{0}$)
has a strongly convex objective function in the decision variables
$(\\{\phi_{i}\\}_{1}^{n},\\{\bm{\xi}_{i}\\}_{1}^{n})$ — its Lagrangian dual
(see ($D$)) is a convex QP with $O(n^{2})$ variables over the nonpositive
orthant. We present large scale algorithms for this dual and study their
computational guarantees. The large number of variables poses computational
challenges for full-gradient-based optimization methods as soon as $n$ becomes
larger than a few thousand. However, we anticipate that $O(n)$-many of the
constraints in ($P$) will suffice. We draw inspiration from the one
dimensional case ($d=1$), an observation that was also used by [6]. Hence, we
use methods inspired by constraint generation [7], which we also refer to
(with an abuse of terminology) as active set type methods [5]. Every step of
our algorithm considers a reduced dual problem where the decision variables,
informally speaking, correspond to a subset of the primal constraints. The
vanilla version of this active set method, which solves the reduced dual
problem to optimality111This is similar to the method of [6] who consider a
constraint generation method for (2) where the reduced sub-problems are solved
to optimality., becomes expensive when $n$ and/or the size of the active set
becomes large especially if one were to perform several active-set iterations.
We propose improved algorithms that perform inexact (or approximate)
optimization for this reduced problem initially and then increase the
optimization accuracy at a later stage. Upon solving the reduced problem
(exactly or inexactly), we examine optimality conditions for the full problem;
and include additional variables into the dual problem, if necessary. To
augment the current active-set, greedy deterministic augmentation rules that
scan all $O(n^{2})$-constraints, become a computational bottleneck —
therefore, we use randomized rules, which leads to important cost savings.
These randomized augmentation rules extend the random-then-greedy selection
strategies proposed by [24] in the context of Gradient Boosting Machines [11].
Our approach operates on the dual and results in a dual feasible solution—we
show how to obtain a primal feasible solution, delivering a duality gap
certificate.

We establish a novel linear convergence rate of our algorithm on the dual,
which is not strongly convex. Our guarantees apply to both exact/inexact
optimization of the reduced problem; and both deterministic and randomized
augmentation rules. As we focus on large scale problems (e.g. $n\geq 10,000$),
inexact optimization of the reduced sub-problem and randomized augmentation
rules play a key role in computational efficiency. As we carefully exploit
problem-structure, our standalone algorithms enjoy a low memory footprint and
can approximately solve convex regression problems with $n\approx 10^{5}$ and
$d=10$ (i.e., QPs with $O(10^{5})$ variables and $O(10^{10})$ constraints) in
minutes. Numerical comparisons suggest that our approach can outperform
existing approaches by a significant margin.

Organization of paper: Section 2 presents both primal and dual formulations of
the full problem ($P_{0}$); and a first order method on the dual. Section 3
presents active-set type algorithms, augmentation rules and associated
computational guarantees. Section 4 discusses computing duality gap
certificates. Section 5 presents numerical experiments. Some technical details
are relegated to Appendix A to improve readability.

Notations: For convenience, we list some notations used throughout the paper.
We denote the set $\\{1,2,\ldots,n\\}$ by $[n]$. The cardinality of a set $W$
is denoted by $|W|$. We denote by $\mathbb{R}_{+}$, $\mathbb{R}_{++}$ the set
of nonnegative and positive real numbers (respectively). A similar notation
applies for $\mathbb{R}_{-}$ and $\mathbb{R}_{--}$. Symbols $\bm{1}_{n}$,
$\bm{e}_{i}$ and $\bm{I}$ denote: a vector of length $n$ of all ones, the
$i$-th standard basis element and the identity matrix (respectively).
$\mathrm{span}(A)$ denotes the linear space generated by the vectors in the
set $A$. For a matrix $\bm{B}$, let $\mathrm{vec}(\bm{B})$ denote a vectorized
version of $\bm{B}$. The largest singular value of a matrix $\bm{B}$ is
denoted by $\lambda_{\max}(\bm{B})$. We use $\lVert\cdot\rVert$ to denote the
Euclidean norm of a vector and the spectral norm of a matrix. A convex
function $f$ is said to be $\sigma$-smooth if it is continuously
differentiable with $\sigma$-Lipschitz gradient; $f$ is $\mu$-strongly convex
if $f(\bm{x})-\frac{\mu}{2}\bm{x}^{\top}\bm{x}$ is convex. Finally, $\partial
f(\bm{x})$ denotes the set of subgradients of $f$ at $\bm{x}$.

##  2 Primal and Dual Formulations

We introduce some notation to write Problem ($P_{0}$) compactly. Let
$\bm{y}=[y_{1},\ldots,y_{n}]^{\top}\in\mathbb{R}^{n}$,
$\bm{X}=[\bm{x}_{1}^{\top},\ldots,\bm{x}_{n}^{\top}]$ $\in\mathbb{R}^{n\times
d}$, $\bm{\phi}=[\phi_{1},\ldots,\phi_{n}]^{\top}\in\mathbb{R}^{n}$, and
$\bm{\xi}=\mathrm{vec}([\bm{\xi}_{1},\ldots,\bm{\xi}_{n}])\in\mathbb{R}^{nd}$.
We define $\bm{A}\in\mathbb{R}^{n(n-1)\times n}$ and
$\bm{B}\in\mathbb{R}^{n(n-1)\times nd}$ such that the rows of
$\bm{A}\bm{\phi}+\bm{B}\bm{\xi}$ correspond to the constraints
$\phi_{j}-\phi_{i}-\langle\bm{x}_{j}-\bm{x}_{i},\bm{\xi}_{i}\rangle$ for
$(i,j)\in\Omega$. Hence ($P_{0}$) is equivalent to:

$\displaystyle\begin{aligned}
&\underset{\bm{\phi},\bm{\xi}}{\text{minimize}}&&f(\bm{\phi},\bm{\xi}):=\frac{1}{2}\lVert\bm{y}-\bm{\phi}\rVert^{2}+\frac{\rho}{2}\lVert\bm{\xi}\rVert^{2}~{}~{}~{}~{}~{}\text{s.t.}&&\bm{A\phi}+\bm{B\xi}\geq\bm{0},\end{aligned}$
($P$)

where, $\bm{A\phi}+\bm{B\xi}\geq\bm{0}$ denotes componentwise inequality. Due
to strong convexity of the objective, ($P$) has a unique minimizer denoted by
$(\bm{\phi}^{\star},\bm{\xi}^{\star})$. While ($P$) has a large number (i.e.,
$O(n^{2})$) of constraints, given that the affine hull of the data points has
full dimension222This occurs with probability one if the covariates are drawn
from a continuous distribution., i.e.
$\mathrm{span}(\\{\bm{x}_{i}-\bm{x}_{j}\\}_{j\neq i})=\mathbb{R}^{d}$, it can
be shown that the constraint matrix
$\bm{C}=\begin{bmatrix}\bm{A}&\bm{B}\end{bmatrix}\in\mathbb{R}^{n(n-1)\times(n+nd)}$
is of rank $O(nd)$. This is a key motivation for our active-set approach.

The Lagrangian dual of ($P$) is equivalent to the following convex problem

$L^{\star}=\underset{\bm{\lambda}\in\mathbb{R}^{n(n-1)}}{\mathrm{minimize}}\quad
L(\bm{\lambda}):=\frac{1}{2\rho}\bm{\lambda}^{\top}(\rho\bm{AA}^{\top}+\bm{BB}^{\top})\bm{\lambda}-\bm{y}^{\top}\bm{A}^{\top}\bm{\lambda},\quad\mathrm{s.t.}\quad\bm{\lambda}\leq
0.$ ($D$)

Note that $\bm{\lambda}\mapsto L(\bm{\lambda})$ is not strongly convex, but it
is $\sigma$-smooth, where $\sigma$ is the maximum eigenvalue of
$\bm{Q}:=\bm{A}\bm{A}^{\top}+\frac{1}{\rho}\bm{B}\bm{B}^{\top}$.

Unlike the primal ($P$), the dual problem ($D$) is amenable to proximal
gradient methods [29, 4] (PGD). Other gradient based methods like accelerated
proximal gradient methods (APG) [4, 29], the limited-memory Broyden-Fletcher-
Goldfarb-Shanno (L-BFGS) method [23] (for example), may also be used; and they
work well in our numerical experience. However, every iteration of PGD
requires computing the gradient $\nabla
L(\bm{\lambda})\in\mathbb{R}^{n(n-1)}$. While an unstructured gradient
computation will cost $O(n^{4})$, exploiting the structure of $\bm{A},\bm{B}$,
this cost can be reduced to $O(n^{2}d)$, allowing us to scale these algorithms
for instances with $n\approx 3,000$. Such matrix-vector multiplications can be
used to estimate $\sigma$ via the power method or backtracking line-search
[4]. Under the Lipschitz continuity of $\nabla L$, the standard sublinear
convergence rate $O(1/t)$ holds [4]. With additional strong convexity
assumption, PGD converges at a linear rate [28]. Although $L(\bm{\lambda})$ is
not strongly convex, it satisfies the Polyak-Łojasiewicz condition [19, 27]
under which it enjoys a linear convergence rate.

##  3 Active Set Type Algorithms

As ($D$) has $O(n^{2})$ variables, the proximal gradient method (owing to full
gradient computations) becomes prohibitively expensive when $n$ becomes larger
than a few thousand. However, as discussed earlier, we expect $O(n)$-many of
these variables to be nonzero at an optimal solution—motivating the use a
constraint-generation/active set-type method on the primal ($P$), which
relates to a column-generation type method on the dual ($D$).

Constraint generation is traditionally used in the context of solving large-
scale linear programs [10, 7]. When used in the context of the QP ($P$), as
done in [6], we start with a reduced problem with a small subset of
constraints in ($P$). With a slight abuse of nomenclature, we refer to these
constraints as the active set333Our use of the term “active set” differs from
the active set method for solving a QP, as discussed in [30]. After obtaining
an optimal dual solution to this reduced problem, the traditional form of
constraint generation will augment the active set with some dual variables
that correspond to the violated primal constraints (if any) and re-solve the
problem on the expanded set of constraints. We mention a couple of
shortcomings of this approach: (a) Solving the reduced problem to optimality
becomes expensive (especially when the active set becomes large and/or if
several iterations of constraint-generation is needed); and (b) finding the
variables that need to be appended to the active set has a large cost of
$O(n^{2}d)$ operations.

To circumvent these shortcomings, we propose modifications to the above
constraint generation or active set method. To address (a), we solve the
reduced sub-problem inexactly (e.g., by taking a few iterations of the
proximal gradient method). To address (b), we consider randomized rules to
reduce the cost of augmenting the active set from $O(n^{2}d)$ to $O(nd)$ (for
example). We show that our proposed algorithm converges; and does so with a
linear convergence rate.

###  3.1 Properties of the reduced problem

Let $W\subseteq\Omega$ index a subset of the constraints in ($P$); and
consider the reduced primal:

$\displaystyle\begin{aligned}
&\underset{\bm{\phi},\bm{\xi}}{\text{minimize}}&&f(\bm{\phi},\bm{\xi})=\frac{1}{2}\lVert\bm{y}-\bm{\phi}\rVert^{2}+\frac{\rho}{2}\lVert\bm{\xi}\rVert^{2}~{}~{}~{}~{}~{}\text{s.t.}&&\bm{A}_{W}\bm{\phi}+\bm{B}_{W}\bm{\xi}\geq\bm{0}\end{aligned}$
($P_{W}$)

where, $\bm{A}_{W}$ (and $\bm{B}_{W}$) denotes matrix $\bm{A}$ (and $\bm{B}$)
restricted to rows indexed by $W$.

In the rest of the paper, we will use $W$ as a subscript for vectors or
matrices whose size changes with $W$, and use $W$ (or $[W]$) as superscript
for vectors or matrices whose size does not change with $W$. When $W=\Omega$,
the relaxed problem is the original problem, and we drop the use of $\Omega$
as subscript and/or superscript for notational convenience.

We consider solving the dual of ($P_{W}$). Proposition 1 presents some of its
properties.

######  Proposition 1.

The Lagrangian dual of ($P_{W}$) is given by:

$\underset{\bm{\lambda}_{W}}{\mathrm{min}}~{}~{}L_{W}(\bm{\lambda}_{W}):=\tfrac{1}{2\rho}\bm{\lambda}_{W}^{\top}(\rho\bm{A}_{W}\bm{A}_{W}^{\top}+\bm{B}_{W}\bm{B}_{W}^{\top})\bm{\lambda}_{W}-\bm{y}^{\top}\bm{A}_{W}^{\top}\bm{\lambda}_{W}~{}~{}\mathrm{s.t.}~{}~{}\bm{\lambda}_{W}\leq\bm{0}.$
($D_{W}$)

Let $L_{W}^{\star}$ be the optimal objective value for ($D_{W}$). The
objective function $L_{W}(\cdot):\mathbb{R}^{|W|}\to\mathbb{R}$ is
$\sigma_{W}$-smooth for some $\sigma_{W}\leq\sigma$. The set of all optimal
solutions to ($D_{W}$) is a polyhedron of the form

$\Lambda^{\star}_{W}=\left\\{\bm{\lambda}_{W}\in\mathbb{R}^{|W|}:\bm{A}_{W}^{\top}\bm{\lambda}_{W}=\bm{s}_{A}^{W},\frac{1}{\sqrt{\rho}}\bm{B}_{W}^{\top}\bm{\lambda}_{W}=\bm{s}_{B}^{W},\bm{\lambda}_{W}\leq
0\right\\}$ (3)

where,
$\bm{s}_{A}^{W}=\bm{y}-\bm{\phi}^{W}_{\star},\bm{s}_{B}=-\sqrt{\rho}\bm{\xi}^{W}_{\star}$;
and $(\bm{\phi}_{\star}^{W},\bm{\xi}_{\star}^{W})$ is the unique optimal
solution to the reduced primal problem ($P_{W}$).

###### Proof Sketch.

The Lagrangian dual is obtained by computing
$\min_{\bm{\xi},\bm{\phi}}\mathcal{L}(\bm{\phi},\bm{\xi};\bm{\lambda}_{W})$
where,
$\mathcal{L}(\bm{\phi},\bm{\xi};\bm{\lambda}_{W})=f(\bm{\phi},\bm{\xi})+\bm{\lambda}_{W}^{\top}(\bm{A}_{W}\bm{\phi}+\bm{B}_{W}\bm{\xi})$
for any $\bm{\lambda}_{W}\leq 0$ and then flipping the sign of the function to
get a convex function. As $\bm{\lambda}_{W}\mapsto L_{W}(\bm{\lambda}_{W})$ is
a quadratic function, it is $\sigma_{W}$-smooth with $\sigma_{W}$ being the
maximum eigenvalue of the Hessian matrix. Since the Hessian of the reduced
dual problem is a submatrix of that of the original dual problem, it follows
that $\sigma_{W}\leq\sigma$. The formulation of the polyhedral set follows
from the KKT conditions [9]. ∎

PGD is readily applicable to ($D_{W}$). The per iteration cost is $O(|W|d)$,
which is much smaller than $O(n^{2}d)$ (as $|W|\sim n$). Solving the reduced
problem ($D_{W}$) is usually much faster than solving the full problem ($D$)
when $|W|\ll|\Omega|$.

###  3.2 Augmentation Rules

For any $W\subseteq\Omega$, given a feasible solution444This can be obtained
by solving ($D_{W}$) exactly (i.e., to optimality) or inexactly (i.e.,
approximately). $\bm{\lambda}_{W}\in\mathbb{R}_{-}^{|W|}$ to the reduced dual
($D_{W}$), we can construct its corresponding primal variables
$(\bm{\phi}^{W},\bm{\xi}^{W})$ for ($P_{W}$), by making use of the KKT
conditions:

$(\bm{\phi}^{W},\bm{\xi}^{W})=(\bm{y}-\bm{A}^{\top}_{W}\bm{\lambda}_{W},-\frac{1}{\rho}\bm{B}_{W}^{\top}\bm{\lambda}_{W})\in\mathbb{R}^{n}\times\mathbb{R}^{nd}.$
(4)

Notice that for a general dual variable
$\bm{\lambda}_{W}\in\mathbb{R}^{|W|}_{-}$, the primal variables obtained via
(4) may not be feasible for the reduced primal problem ($P_{W}$). Below we
discuss some rules for augmenting the current set of constraints (i.e., the
active set).

After obtaining $(\bm{\phi}^{W},\bm{\xi}^{W})$, we check if it is feasible for
($P$)—that is, we verify if each component of $\bm{A\phi}^{W}+\bm{B\xi}^{W}$
(denoted by $v_{(i,j)}$) is nonnegative:

$v_{(i,j)}=\phi_{j}^{W}-\phi_{i}^{W}-\langle\bm{x}_{j}-\bm{x}_{i},\bm{\xi}_{i}^{W}\rangle\geq
0~{}~{}~{}\forall(i,j)\in\Omega.$ (5)

To this end, it is helpful to decompose the $O(n^{2})$ constraints into $n$
blocks

$\Omega=\bigcup_{i}\Omega_{i\cdot}=\bigcup_{j}\Omega_{\cdot j},$

where $\Omega_{i\cdot}=\\{(i,j):j\neq i,1\leq j\leq n\\}$ and $\Omega_{\cdot
j}=\\{(i,j):i\neq j,1\leq i\leq n\\}.$ Similarly, we define the slices
$W_{i\cdot}$ and $W_{\cdot i}$ for all $i$. The two decompositions
$\\{\Omega_{i\cdot}\\}_{1}^{n}$ and $\\{\Omega_{\cdot j}\\}_{1}^{n}$ have
different geometric interpretations. Now define points
$P_{j}=\\{\bm{x}_{j},\phi_{j}^{W}\\}$ and hyperplanes
$H_{i}:y=\langle\bm{x}-\bm{x}_{i},\bm{\xi}^{W}_{i}\rangle+\phi_{i}^{W}$ in
$\mathbb{R}^{d+1}$. Note that each point $P_{i}$ lies on the hyperplane
$H_{i}$, and $v_{(i,j)}$ denotes the vertical distance between $P_{j}$ and
$H_{i}$. For each $i$, the nonnegativity of $v_{(i,j)}$ for all
$(i,j)\in\Omega_{i\cdot}$, checks if the hyperplane $H_{i}$ lies below each
point $P_{j}$, i.e. if $H_{i}$ is a supporting hyperplane of the points
$\\{P_{j}\\}_{1}^{n}$. On the other hand, for each $j$, the nonnegativity of
$v_{(i,j)}$ for all $(i,j)\in\Omega_{\cdot j}$, checks if $P_{j}$ lies above
each hyperplane $H_{i}$—i.e., if $P_{j}$ lies above the piecewise maximum of
these hyperplanes. The quantity $|\min_{i}v_{(i,j)}|$ is the amount by which
$P_{j}$ lies below the piecewise maximum.

In Section 4 we see that the block decomposition interpretation is useful in
obtaining a primal feasible solution. Next we discuss a deterministic
augmentation rule—due to its high computational cost, we subsequently present
randomized augmentation rules—both of which make use of the above
decomposition of $\Omega$.

####  3.2.1 Deterministic Augmentation Rule

We first present a simple greedy-like deterministic augmentation rule:

  * Rule 1.

Greedy within each Block: For each block $\Omega_{i\cdot}$ (or $\Omega_{\cdot
j}$), choose $P$ pairs with the smallest $v_{(i,j)}$-values among
$\Omega_{i\cdot}\backslash W_{i\cdot}$ (or $\Omega_{\cdot j}\backslash
W_{\cdot j}$). From these $P$ indices, we add to the current active set $W$,
only those with negative $v_{(i,j)}$-values.

Rule 1 evaluates $O(n^{2})$ constraints with computational cost
$O(n^{2}d+d|W|)$ and augments $W$ by at most $nP$ such constraints. For every
block, obtaining the largest $P$ violations can be done via a partial
sort—leading to a total cost of $O(n^{2}d+nP\log P)$ for $n$ blocks555A
similar augmentation rule with $P=1$ is used in [6]. As a result, with this
greedy rule, the augmentation becomes a major computational bottleneck for the
active-set type method. This motivates the randomized augmentation rules
discussed below.

####  3.2.2 Randomized Augmentation Rules

We present four randomized rules that sample a small subset of the indices in
$\Omega\setminus W$ instead of performing a (computationally intensive) full
scan across $O(n^{2})$ indices as in Rule 1.

  * Rule 2.

Random: Uniformly sample $K$ indices in $\Omega\backslash W$—the cost of
computing the corresponding $v_{(i,j)}$-values is $O(Kd+d|W|)$.

  * Rule 3.

Random within each Block: For each block $\Omega_{i\cdot}$ (or $\Omega_{\cdot
j}$), uniformly sample $P$ elements in $\Omega_{i\cdot}\backslash W_{i\cdot}$
(or $\Omega_{\cdot j}\backslash W_{\cdot j}$); the total computational cost is
$O(nPd+d|W|)$.

  * Rule 4.

Random then Greedy: Uniformly sample $M$-many $(i,j)$-pairs from
$\Omega\backslash W$, and from these pairs choose the $K$ pairs with the
smallest $v_{(i,j)}$-values. Computing the $M$ values of $v_{(i,j)}$ cost
$O(Md+d|W|)$, and the greedy selection step costs $M+K\log K$. The total
computational cost is $O(Md+d|W|)$.

  * Rule 5.

Random Blocks then Greedy within each Block: Uniformly sample $G$ groups from
$\\{\Omega_{i\cdot}\\}_{i=1}^{n}$ (or $\\{\Omega_{\cdot j}\\}_{j=1}^{n}$) and
for each group, choose the $P$ pairs that have the smallest $v_{(i,j)}$ values
with $(i,j)\in\Omega_{i\cdot}\backslash W_{i\cdot}$ (or $\Omega_{\cdot
j}\backslash W_{\cdot j}$). Similar to Rule 4, the total computational cost is
$O(Gnd+d|W|)$.

Note that the above rules lead to a set of indices denoted by
$\Delta^{\prime}$. From these candidates, we only select those with negative
$v_{(i,j)}$-values, which are then appended to the current active set $W$.
That is, if $\Delta=\\{(i,j)\in\Delta^{\prime}:v_{(i,j)}<0\\}$ then we set
$W\leftarrow W\cup\Delta$.

Of the above, Rules 4 and 5 are inspired by, but different from the random-
then-greedy coordinate descent procedure [24], proposed in the context of
Gradient Boosting.

####  3.2.3 Norms Associated with the Augmentation Rules

The computational guarantees of our active set-type algorithms depend upon
certain norms induced by the aforementioned augmentation rules. We first
present some notation that we will use in our analysis.

######  Definition 1.

Given a vector $\bm{\theta}$, an index set $S$, and $k\leq|S|$, let
$\mathcal{G}[S,k,\bm{\theta}]$ be the set of $k$ elements with the largest
values of $|\theta_{\omega}|$ for $\omega\in S$. Given a pair $(S,k)$ as
above, we let $\mathcal{U}[S,k]$ denote the set of $k$ uniformly subsampled
indices from the set $S$.

######  Definition 2.

Given a vector $\bm{\theta}\in\mathbb{R}^{n(n-1)}$ indexed by $\Omega$, let
$\\{S_{i}\\}_{1}^{n}$ be disjoint subsets of $\Omega$, with
$\bar{S}=\cup_{i\in[n]}S_{i}$. Given positive integers $P,K,M,G$, we define
$\delta_{1}(\bm{\theta},\\{S_{i}\\}),\ldots,\delta_{5}(\bm{\theta},\\{S_{i}\\})\subset\Omega$
as follows:

$\displaystyle\delta_{1}(\bm{\theta},\\{S_{i}\\})=\bigcup\limits_{i=1}^{n}\mathcal{G}[S_{i},P,\bm{\theta}],~{}~{}~{}\delta_{2}(\bm{\theta},\\{S_{i}\\})=\mathcal{U}[\bar{S},K],~{}~{}~{}~{}\delta_{3}(\bm{\theta},\\{S_{i}\\})=\bigcup_{i=1}^{n}\mathcal{U}[S_{i},P]$
$\displaystyle\delta_{4}(\bm{\theta},\\{S_{i}\\})=\mathcal{G}[\mathcal{U}[\bar{S},M],K,\bm{\theta}],~{}~{}~{}~{}~{}~{}~{}\delta_{5}(\bm{\theta},\\{S_{i}\\})=\bigcup_{i\in\mathcal{U}[[n],G]}\mathcal{G}[S_{i},P,\bm{\theta}].$

With the above notation, we can express the indices to be augmented as a
function of the violations ${\bm{v}}$, a vectorized representation of the
entries $\\{v_{(i,j)}\\}_{(i,j)\in\Omega}$ (5). Let
$\Delta^{\prime}_{\\{\ell\\}}$ denote the pairs selected by Rule $\ell$, and
$\Delta_{\\{\ell\\}}=\\{\omega\in\Delta_{\\{\ell\\}}^{\prime}:v_{(i,j)}<0\\}\subseteq\Delta_{\\{\ell\\}}^{\prime}$
be the final set to be added to $W$. Let $\tilde{\bm{v}}$ be a vector having
the same length as $\bm{v}$, with its entries given by
$\tilde{v}_{(i,j)}=\min\\{v_{(i,j)},0\\}$. Then it is easy to verify that
$\Delta^{\prime}_{\\{\ell\\}}$ can be written as
$\delta_{\ell}(\tilde{\bm{v}},\\{S_{i}\\})$ with
$S_{i}=\Omega_{i\cdot}\setminus W_{i\cdot}$ or $S_{i}=\Omega_{\cdot
i}\setminus W_{\cdot i}$.

######  Definition 3.

Given $\bm{\theta}\in\mathbb{R}^{n(n-1)}$, let $\\{\Omega_{i}\\}_{1}^{n}$ be
either $\\{\Omega_{i\cdot}\\}_{1}^{n}$ or $\\{\Omega_{\cdot i}\\}_{1}^{n}$.
Define

$\lVert\bm{\theta}\rVert_{\\{\ell\\}}=\left(\mathbb{E}\left[\sum_{\omega\in\delta_{\ell}(\bm{\theta},\\{\Omega_{i}\\})}|\theta_{\omega}|^{2}\right]\right)^{1/2}~{}~{}~{}\text{for}~{}~{}~{}\ell\in\\{1,\ldots,5\\}.$
(6)

For $\ell=1$, there is no randomness in $\delta_{\ell}$, so the expectation
can be removed. For $\ell\geq 2$, the expectation is taken over the randomness
arising from the selection operator $\mathcal{U}$ (cf Definition 1).

Lemma 1 shows that the expression in display (6) leads to a norm on
$\mathbb{R}^{n(n-1)}$. Furthermore, the norm equivalence constants appearing
in (7) determine the convergence rates of Algorithm 1 (see Theorem 1).

######  Lemma 1.

$\lVert\bm{\theta}\rVert_{\\{\ell\\}}$ defined in (6) is a norm on
$\mathbb{R}^{n(n-1)}$. Furthermore, there are constants
$\alpha_{\\{\ell\\}}\geq\beta_{\\{\ell\\}}\geq 1$, listed in Table 1, such
that for all $\bm{x}\in\mathbb{R}^{n(n-1)}$

$\alpha_{\\{\ell\\}}\lVert\bm{x}\rVert_{\\{\ell\\}}^{2}\geq\lVert\bm{x}\rVert_{2}^{2}\geq\beta_{\\{\ell\\}}\lVert\bm{x}\rVert_{\\{\ell\\}}^{2}.$
(7)

The proof of Lemma 1 appears in Section A.1.

Table 1 presents a summary of some key characteristics of the five rules.

Rule | $\Delta^{\prime}_{\\{\ell\\}}$ | $|\Delta^{\prime}_{\\{\ell\\}}|$ | Augmentation Cost | $\alpha_{\\{\ell\\}}$ | $\beta_{\\{\ell\\}}$
---|---|---|---|---|---
1 | $\cup_{i=1}^{n}\mathcal{G}[\Omega_{i},P,\tilde{\bm{v}}]$ | $nP$ | $O(n^{2}d+d|W|)$ | $\frac{n-1}{P}$ | 1
2 | $\mathcal{U}[\Omega,K]$ | $K$ | $O(Kd+d|W|)$ | $\frac{n(n-1)}{K}$ | $\frac{n(n-1)}{K}$
3 | $\cup_{i=1}^{n}\mathcal{U}[\Omega_{i},P]$ | $nP$ | $O(nPd+d|W|)$ | $\frac{n-1}{P}$ | $\frac{n-1}{P}$
4 | $\mathcal{G}[\mathcal{U}[\Omega,M],K,\tilde{\bm{v}}]$ | $K$ | $O(Md+d|W|)$ | $\frac{n(n-1)}{K}$ | $\frac{n(n-1)}{M}$
5 | $\cup_{i\in\mathcal{U}[[n],G]}\mathcal{G}[\Omega_{i},P,\tilde{\bm{v}}]$ | $GP$ | $O(Gnd+d|W|)$ | $\frac{n(n-1)}{GP}$ | $\frac{n}{G}$
Table 1: Summary of some properties of the augmentation rules. Recall that
$\Delta^{\prime}_{\\{\ell\\}}$ denotes the pairs selected as per Rule $\ell$.
The number of candidates to be augmented to the current active set depends
upon the signs of $v_{(i,j)}$s; and is of size at most
$|\Delta^{\prime}_{\\{\ell\\}}|$. Here, Augmentation Cost is the cost of
obtaining $\Delta^{\prime}_{\\{\ell\\}}$. We present estimates of the norm-
equivalence constants $\alpha_{\\{\ell\\}},\beta_{\\{\ell\\}}$ (7). For
notational convenience, we assume that $W=\varnothing$—for a nonempty $W$, we
can replace $\Omega_{i}$ with $\Omega_{i}\backslash W_{i}$.

###  3.3 Active set method with inexact optimization of sub-problems

Once we augment the active set (based on Rules 1–5), we solve the updated
reduced QP either exactly or inexactly. We then identify additional
constraints to be added to the current active set; and continue till some
convergence criteria is satisfied. Our algorithm is summarised below:

Algorithm 1: An Active Set Type Method on the Dual ($D$) Initialize with
$W^{0}$ and $\bm{\lambda}^{0}$. Perform the following steps 1–3 till
convergence. 1. Obtain a feasible solution $\bm{\lambda}_{W^{m}}$ for
$(D_{W^{m}})$ by either (a) solving $(D_{W^{m}})$ to optimality (i.e.,
exactly) or (b) solving $(D_{W^{m}})$ inexactly via a few steps of proximal
gradient descent666Other choices are also possible, as discussed in Remark 2..
2. Compute
$\bm{\phi}^{W^{m}}=\bm{y}-\bm{A}^{\top}_{W^{m}}\bm{\lambda}_{W^{m}},$
$\bm{\xi}^{W^{m}}=-\frac{1}{\rho}\bm{B}_{W^{m}}^{\top}\bm{\lambda}_{W^{m}}$ as
per (4). 3. Use one of the Rules 1–5 to augment $W^{m}$ to obtain $W^{m+1}$;
and go to Step 1.

In what follows, we call Algorithm 1 with option (a) [Step 1] the Exact Active
Set method (EAS); and option (b) [Step 1] the Active Set Gradient Descent
(ASGD) method.

If the size of the active set remains sufficiently small, interior point
solvers (including heavy-duty commercial solvers like Gurobi, Mosek) may be
used for full minimization of the reduced problem, as long as the number of
active set iterations remain small. However, for larger problems (e.g., when
$n\approx 10^{4}$–$10^{5}$), first order methods may be preferable due to
warm-start capabilities (across active sets) and low memory requirements (by
exploiting problem-structure). Furthermore, they also allow for inexact
computation for the sub-problems—see numerical experiments in Section 5
showing the important benefits in approximate optimization.

Theorem 1 establishes that Algorithm 1 requires at most
$m=O(\log(1/\epsilon))$-many outer iterations to deliver an $\epsilon$-optimal
dual solution for Problem ($D$):

######  Theorem 1.

For a given $\bm{\lambda}^{0}$ and $W^{0}$, let
$\bm{\lambda}^{m}=\bm{\lambda}^{[W_{m}]}$ (for $m\geq 1$) be the sequence
generated by Algorithm 1. Then, for either EAS or ASGD with augmentation Rule
$\ell$, we have:

$\mathbb{E}[L(\bm{\lambda}^{m})-L^{\star}]\leq\left(1-\frac{\mu}{\sigma\alpha_{\\{\ell\\}}}\right)^{m}[L(\bm{\lambda}^{0})-L^{\star}],$

where $\sigma$ is the smoothness constant of $L$, $\mu$ is a constant
depending only on $\bm{A},\bm{B}$ and $\rho$, and $\alpha_{\\{\ell\\}}$
appears in Lemma 1.

The proof of Theorem 1 is provided in Sections A.2 and A.3. We make a few
remarks regarding Theorem 1.

######  Remark 1.

If we used PGD on the dual ($D$), the linear convergence rate parameter [29,
19] would be $(1-\mu/\sigma)$. The parameter in Theorem 1 has an additional
factor $1/\alpha_{\\{\ell\\}}$, which arises due to the use of an active set
method in conjunction with the augmentation rules. We present an example in
Section B.1 showing that this factor is unavoidable.

######  Remark 2.

Algorithm 1 allows for both exact and inexact optimization of the reduced
problem—the guarantees for Theorem 1 apply to both variants. Furthermore,
Theorem 1 applies to both deterministic and randomized augmentation rules.
Although our theory is based on a proximal gradient update on the reduced
problem, the theory also applies to other descent algorithms that have a
sufficient decrease condition, e.g. accelerated proximal gradient (APG)
methods [4, 29] and L-BFGS [23] with some modifications–see discussions in
Remark 4.

######  Remark 3.

Theorem 1 implies that we need at least
$O\left(({\sigma\alpha_{\\{\ell\\}}}/{\mu})\log({1}/{\epsilon})\right)$-many
outer iterations, to achieve an accuracy level of $\epsilon>0$. This
quantifies the worst-case convergence rate via $\alpha_{\\{\ell\\}}$. By
following the arguments in the proof, one can obtain a more optimistic upper
bound on the number of iterations, by replacing $\alpha_{\\{\ell\\}}$ with
$\beta_{\\{\ell\\}}$ (cf Lemma 1), as we discuss below.

Understanding the costs of different rules: To better understand the
computational costs of the different rules, we consider a setting where the
maximum number of selected pairs (i.e., $|\Delta_{\\{\ell\\}}^{\prime}|$) is
set to be $O(n)$ across all rules777The choice of $P,K,M,G$ are specified in
the caption of Table 2.

According to Remark 3, to achieve $\epsilon$ accuracy, we need
$O((\sigma\alpha_{\\{\ell\\}}/{\mu})\log({1}/{\epsilon}))$ outer iterations in
the worst case and $O(({\sigma\beta_{\\{\ell\\}}}/{\mu})\log({1}/{\epsilon}))$
outer iterations in the best case. Hence, for a fixed $\epsilon>0$, the total
augmentation cost is proportional to $\alpha_{\\{\ell\\}}\times$Augmentation
Cost for the worst case and $\beta_{\\{\ell\\}}\times$Augmentation Cost for
the best case. See Table 2 for an illustration.

Rule | $|\Delta^{\prime}_{\\{\ell\\}}|$ | Augmentation Cost | $\alpha_{\\{\ell\\}}$ | $\beta_{\\{\ell\\}}$
---|---|---|---|---
1 | $O(n)$ | $O(n^{2}d+d|W|)$ | $O(n)$ | $O(1)$
2 | $O(n)$ | $O(nd+d|W|)$ | $O(n)$ | $O(n)$
3 | $O(n)$ | $O(nd+d|W|)$ | $O(n)$ | $O(n)$
4 | $O(n)$ | $O(n\sqrt{n}d+d|W|)$ | $O(n)$ | $O(\sqrt{n})$
5 | $O(n)$ | $O(n\sqrt{n}d+d|W|)$ | $O(n)$ | $O(\sqrt{n})$

Table 2: Comparing Augmentation Rules: We present an instance of Table 1 where
$|\Delta_{\\{\ell\\}}^{\prime}|$ is the same across $\ell$. Specifically, we
set $P=1$ for Rule 1, $K=n$ for Rule 2, $P=1$ for Rule 3, $M=n\sqrt{n},K=n$
for Rule 4, and $G=\sqrt{n},P=\sqrt{n}$ for Rule 5. Here, we ignore log-terms
in the big $O$ notation.

We make a note of the following observations from Table 2.

  * •

As the maximum size of the augmentation set (i.e,
$|\Delta_{\\{\ell\\}}^{\prime}|$) at each iteration is the same across all
rules, $\alpha_{\\{\ell\\}}$ is the same across all the five rules $\ell$.

  * •

In the worst case, different rules take the same number of iterations to
achieve $\epsilon$ accuracy. The largest total augmentation costs (i.e.,
$\alpha_{\\{\ell\\}}\times$Augmentation Cost) are incurred by the purely
greedy rule (Rule 1), and then the random-then-greedy rules (Rule 4 and Rule
5).

  * •

In the best case, the greedy rules take fewer iterations—the purely greedy
rule (Rule 1) is better than random-then-greedy rules (Rules 4 and 5). The
purely random rules (Rules 2 and 3) have similar iteration counts for the best
and worst cases (as $\alpha_{\\{\ell\\}}=\beta_{\\{\ell\\}}$).

The above discussion pertains to our theoretical guarantees. Our numerical
experiments offer a refined understanding of the operating characteristics of
the different algorithms.

####  3.3.1 Related work

As discussed earlier, Algorithm 1 (with exact optimization for the
subproblems) is inspired by constraint generation, a classic tool for solving
large scale linear programs [10, 7]. In the context of convex regression,
constraint generation ideas are used in [6] but our framework differs. [6] use
a primal approach and solve the restricted primal problem to optimality using
commercial QP solvers (e.g., Gurobi). They do not present computational
guarantees for their procedure. In this paper we consider the subgradient
regularized problem ($P$), and our algorithms are on the dual problem ($D$).
We extend the framework of [6] to allow for both exact and inexact
optimization of the reduced problems. As a main computational bottleneck in
the constraint generation procedure is in the augmentation step which requires
scanning all the $O(n^{2})$ constraints, we present randomized augmentation
rules, which examine only a small subset of the constraints, before augmenting
the active-set. Our computational guarantees, which to our knowledge are
novel, apply to both the exact and inexact sub-problem optimization processes
and deterministic/randomized augmentation rules.

As far as we can tell, with the exception of [6], convex optimization-based
approaches for convex regression (e.g, [25, 22, 1]) that precede our work, do
not use active set-type methods and hence apply to smaller problem instances
$n\approx 3000$.

We emphasize that the main focus of our paper is to approximately solve large
scale instances of Problem ($D$) (e.g., $n\approx 10^{4}$ to $10^{5}$), and
numerical experiments presented in this paper suggest that our proposal has an
edge over competing approaches.

##  4 Primal Feasibility and Duality Gap

Given a feasible solution for the full dual ($D$), we show how to obtain a
feasible solution for the primal ($P$) (note that (4) need not deliver a
primal feasible solution). A primal feasible solution is useful as it leads to
a convex function estimate and also a duality gap certificate.

Given a primal candidate $(\bm{\phi},\bm{\xi})$, we will construct
$(\tilde{\bm{\phi}},\tilde{\bm{\xi}})$ that is feasible for ($P$). To this
end, let
$v_{(i,j)}=\phi_{j}-\phi_{i}-\langle\bm{x}_{j}-\bm{x}_{i},\bm{\xi}_{i}\rangle,$
for all $i\neq j$ and we use the convention $v_{(j,j)}=0$. Furthermore, let

$\nu_{j}=\min_{i}~{}v_{(i,j)}~{}~{}~{}\text{and}~{}~{}~{}\kappa_{j}\in\arg\min\\{\lVert\bm{\xi}_{k}\rVert:k\in\arg\min_{i}v_{(i,j)}\\}$
(8)

and define
$(\tilde{\bm{\phi}},\tilde{\bm{\xi}}):=\bm{\psi}(\bm{\phi},\bm{\xi})$ as
follows

$\tilde{\bm{\xi}}_{j}=\bm{\xi}_{\kappa_{j}}~{}~{}~{}\text{and}~{}~{}~{}~{}\tilde{\bm{\phi}}=\bm{\phi}-\bm{\nu}+c\bm{1}$
(9)

where, $c$ is such that $\bm{1}^{\top}\tilde{\bm{\phi}}=\bm{1}^{\top}\bm{y}$.
The following proposition shows that $(\tilde{\bm{\phi}},\tilde{\bm{\xi}})$ is
feasible for the full primal problem.

######  Proposition 2.

For any pair $(\bm{\phi},\bm{\xi})$, the pair
$(\tilde{\bm{\phi}},\tilde{\bm{\xi}})$ defined in (9), is feasible for ($P$).
Furthermore, if $(\bm{\phi}^{\star},\bm{\xi}^{\star})$ is an optimal solution
for ($P$) then it is a fixed point for the map $\bm{\psi}(\cdot,\cdot)$, i.e.,
$\bm{\psi}(\bm{\phi}^{\star},\bm{\xi}^{\star})=(\bm{\phi}^{\star},\bm{\xi}^{\star}).$

###### Proof of Proposition 2.

Given $\\{{\phi}_{i}\\}$, $\\{\bm{\xi}_{i}\\}$ and a scalar $c$, we define the
following piecewise linear convex function (in $\bm{x}$):

$\varphi(\bm{x})=\max_{i\in[n]}~{}\\{\langle\bm{x},\bm{\xi}_{i}\rangle+(\phi_{i}-\langle\bm{x}_{i},\bm{\xi}_{i}\rangle)\\}+c.$

For any $j$, we have

$\max_{i\in[n]}\\{\langle\bm{x}_{j},\bm{\xi}_{i}\rangle+(\phi_{i}-\langle\bm{x}_{i},\bm{\xi}_{i}\rangle)\\}=\max_{i\in[n]}\\{\phi_{j}-(\phi_{j}-\phi_{i}-\langle\bm{x}_{j}-\bm{x}_{i},\bm{\xi}_{i}\rangle)\\}=\phi_{j}-\min_{i\in[n]}v_{(i,j)}.$

By definition of $\tilde{\bm{\phi}},\tilde{\bm{\xi}}$ given by (9), it follows
that

$\varphi(\bm{x}_{j})=\tilde{\phi}_{j}\quad\text{and}\quad\tilde{\bm{\xi}}_{j}=\bm{\xi}_{\kappa_{j}}\in\partial\varphi(\bm{x}_{j}).$

Therefore, $(\tilde{\bm{\phi}},\tilde{\bm{\xi}})$ is feasible for the full
primal problem.

Due to the feasibility of $(\bm{\phi}^{\star},\bm{\xi}^{\star})$, we have
$\bm{v}^{\star}=\bm{A}\bm{\phi}^{\star}+\bm{B}\bm{\xi}^{\star}\geq\bm{0}$, and
thus $\min_{i}v_{(i,j)}=0=v_{(j,j)}$. This means $\bm{\nu}=\bm{0}$, so it
suffices to show that $c=0$. Note that the first term in the primal objective
$\frac{1}{2}\lVert\bm{\phi}-\bm{y}\rVert^{2}$ is minimized when $\bm{\phi}$
and $\bm{y}$ have the same mean (otherwise we can add a constant to
$\bm{\phi}$ to decrease the objective). Therefore, the $c$ that satisfies
$\bm{1}^{\top}(\bm{\phi}^{\star}+c\bm{1})=\bm{1}^{\top}\bm{y}$ must be 0.
Thus, we have shown that $\tilde{\bm{\phi}}^{\star}=\bm{\phi}^{\star}$. If
$\tilde{\bm{\xi}}^{\star}\neq\bm{\xi}^{\star}$, then there exists $j$ such
that $\kappa_{j}\neq j$, which means
$\lVert\bm{\xi}_{\kappa_{j}}^{\star}\rVert\leq\lVert\bm{\xi}_{j}^{\star}\rVert$,
and we have $f(\tilde{\bm{\phi}}^{\star},\tilde{\bm{\xi}}^{\star})\leq
f(\bm{\phi}^{\star},\bm{\xi}^{\star})$. This contradicts the strict minimality
of $(\bm{\phi}^{\star},\bm{\xi}^{\star})$ as the primal objective function $f$
is strictly convex. ∎

As discussed in Section 3.2, each $(\phi_{i},\bm{\xi}_{i})$ defines a
hyperplane $H_{i}:y=\langle\bm{x}-\bm{x}_{i},\bm{\xi}_{i}\rangle+\phi_{i}$
containing $P_{i}=(\bm{x}_{i},\phi_{i})$. However, $(\bm{\phi},\bm{\xi})$ may
not satisfy the convexity constraint (i.e., $P_{j}$ may not lie above
$H_{i}$); and $|\nu_{j}|$ quantifies how far $P_{j}$ lies below the piecewise
maximum of $H_{i}$s. Intuitively, $(\tilde{\bm{\phi}},\tilde{\bm{\xi}})$
attains feasibility by taking a piecewise maximum of these hyperplanes and
then adjusts itself by a constant in an attempt to decrease the primal
objective.

Duality Gap: In light of strong duality between ($P$) and ($D$), we have
$L(\bm{\lambda}^{\star})=-f(\bm{\phi}^{\star},\bm{\xi}^{\star})$. From a dual
feasible solution $\bm{\lambda}$, we can obtain a primal variable
$(\bm{\phi},\bm{\xi})$ via (4). If $\bm{\psi}(\bm{\phi},\bm{\xi})$ is a primal
feasible solution obtained from $(\bm{\phi},\bm{\xi})$ by the operation
defined in (9), then we can compute a duality gap via:

$\underline{L}(\bm{\lambda})=-f(\bm{\psi}(\bm{\phi},\bm{\xi}))\leq
L^{\star}\leq L(\bm{\lambda}).$

The gap $L(\bm{\lambda})-\underline{L}(\bm{\lambda})$ equals zero if and only
if $\bm{\lambda}$ is an optimal solution for ($D$).

##  5 Numerical Experiments

We present numerical experiments to study the different variants of our
algorithm; and compare it with current approaches. As our focus is on solving
large-scale instances of ($P_{0}$), we present instances (both real and
synthetic datasets) with $n\sim 10^{4}$-$10^{5}$ and $20\geq d\geq 4$.

Datasets: We consider the following synthetic and real datasets for the
experiments.

Synthetic Data: Following [25], we generate data via the model
$y_{i}=\phi^{0}(\bm{x}_{i})+\epsilon_{i},i\in[n]$ where, $\phi^{0}(\bm{x})$ is
a convex function,
$\epsilon_{i}\stackrel{{\scriptstyle\text{iid}}}{{{\sim}}}N(0,\gamma^{2})$ and
$\gamma$ is chosen to match a specified value of signal-to-noise ratio,
$\mathrm{SNR}=\|\bm{\phi}^{0}\|^{2}/\|\bm{\epsilon}\|^{2}=3$. The covariates
are drawn independently from a uniform distribution on $[-1,1]^{d}$. Every
feature is normalized to have zero mean and unit $\ell_{2}$-norm; and $\bm{y}$
has zero mean and unit $\ell_{2}$-norm. We consider the following choices of
the underlying true convex function $\phi^{0}(\bm{x})$.

  * •

SD1: $\phi^{0}(\bm{x})=\lVert\bm{x}\rVert_{2}^{2}$

  * •

SD2: $\phi^{0}(\bm{x})=\max_{1\leq i\leq 2d}\\{\bm{\xi}_{i}^{\top}\bm{x}\\}$,
where $\bm{\xi}_{i}\in\mathbb{R}^{d}$ are independently drawn from the uniform
distribution $[-1,1]^{d}$.

Real Data: We also consider the following real datasets in our experiments:

  * •

RD1: This dataset has been used in earlier work [26, 25, 6]. Here, $n=30,000$
and $d=4$.

  * •

RD2: This dataset taken from [32] and studied in [16, 3], has $n=15,000$ and
$d=4$.

Additional discussions on the real datasets are available in Section C.1. In
both examples above, covariates and response are centered and scaled so that
each variable has unit $\lVert\cdot\rVert_{2}$-norm.

Algorithms: As our main focus is on large scale instances (e.g., $n\geq
10,000$)—we compare our approach (with different rules) versus the state-of-
the-art cutting plane (CP) method [6] using the authors’ implementation. We
also consider Rule 1 with exact optimization of sub-problems (denoted by Rule1
in Figure 1 and Table 3)—this is essentially the same as CP, but with our own
implementation of an AGD method with restarts [31] (we also experimented with
LBFGS, and the runtimes were very similar) in place of Gurobi’s QP solver. In
addition, we also consider a two-stage algorithm with Rules 1, 2 and 4
(denoted by 2Stage-Rule#) where in the first stage we perform inexact
optimization over the active set followed by exact optimization in the second
stage. Empirically, we have observed that inexact optimization methods are
useful in the initial stages, but then the progress becomes slow, as the
augmentation rules add only few constraints in each iteration. Therefore, when
the number of added constraints is less than $0.005n$ for consecutive 5
iterations, we switch to exact optimization888This is a choice we used in our
experiments, and can be tuned in general for performance benefits.. Note that
our theory (Section 3) guarantees convergence for this two-stage procedure.
Exact optimization in the second stage is considerably less expensive compared
to exact optimization in the first stage—this is because, solutions from the
first-stage QPs deliver good warm-starts for the second-stage QPs (the
proximal gradient methods we use for the QPs can use warm-starts). During the
augmentation stage, we consider a constraint to violate feasibility if
$v_{(i,j)}<-\tau$ for some small threshold $\tau=10^{-4}$. For the two-stage
algorithms, we set a tighter violation threshold $\tau=10^{-8}$ in the second
stage. Additional discussions on algorithm parameter choices can be found in
Section C.2.

Software Specifications: All computations were carried out on MIT’s Engaging
Cluster on an Intel Xeon 2.30GHz machine, with one CPU and 8GB of RAM. All of
our algorithms were written in the Julia programming language (v1.5). Our code
is available at

https://github.com/wenyuC94/ConvexRegression

Performance of proposed algorithms:  Figure 1 presents the relative objective
of the dual, defined as

$\text{Rel.
Obj.}={(L(\bm{\lambda}^{t})-L^{\star})}/{\max\\{|L^{\star}|,1\\}},$

for different algorithms versus time (in seconds) and number of active set
iterations. Above, $L^{\star}$ is taken as the minimum objective among
solutions obtained by all the algorithms running for 4 hours. In Figure 1, we
consider synthetic datasets SD1 and SD2 with different $n$’s and $d$’s, as
well as the real dataset RD1. Note that the choices of $\rho$ displayed in the
figure lead to good out-of-sample prediction performance—this is discussed in
further detail in Figure 2. Figure 1 shows that our proposed two-stage active
set algorithm with Rule 4 (i.e., random then greedy rule) performs the best in
both cases. Recall that CP [6] and Rule 1, are essentially the same, except
for different implementations used to solve the reduced QP—we use accelerated
first-order methods, while [6] uses the Gurobi to solve the subproblem. We see
Gurobi runs faster than our algorithm. However, due to the long search time
during the augmentation process, both these methods make little progress in 4
hours.

Table 3 presents a snapshot of the runtimes of our proposed algorithms and CP.
Here, the runtime corresponds to the time (s) taken by an algorithm to achieve
a $0.05$ relative objective. We can see that our proposed algorithm with Rule
2 (random rule) and Rule 4 (random-then-greedy rule) can achieve such accuracy
in 10 minutes, while CP and Rule 1 are unable to achieve such accuracy within
4 hours.

Objective Profiles for SD1, $n=10,000$, $d=4$ and $\rho=10^{-4}$
---
Rel. Obj. |  |
Objective Profiles for SD1, $n=100,000$, $d=20$ and $\rho=10^{-4}$
Rel. Obj. |  |
Objective Profiles for SD2, $n=100,000$, $d=20$ and $\rho=10^{-4}$
Rel. Obj. |  |
Objective Profiles for RD1, $n=30,000$, $d=4$ and $\rho=10^{-5}$
Rel. Obj. |  |

Figure 1: Plots (in log-scale) of Relative Objective versus time (secs) [left
panel]; and # active set iterations [right panel]. We consider three synthetic
data sets (SD1, SD2) with different $(n,d)$-values and $\rho=10^{-4}$, and a
real data set (RD1) with $n=30,000$, $d=4$, and $\rho=10^{-5}$. We consider
different algorithms based on different rules, and compare them against the
cutting plane method [6]. For each algorithm, we run 5 repetitions, each bold
line corresponds to the median of the profiles of one algorithm and each thin
line corresponds to the profile of one repetition. For the left panel, we show
the profiles up to the max time of 4 hours. For the right panel, we show the
profiles up to 250 active set iterations.  Table 3: Comparison of runtime (s)
of our proposed algorithms versus the cutting plane (CP) method of [6]. Here
runtime refers to the time taken to achieve a Rel. Obj. of $5\times 10^{-2}$.
We run 5 replications of each algorithm with different random seeds and report
the median runtime and median absolution deviation (MAD) of runtimes to
achieve the relative objective. ‘*’ means that the algorithm does not achieve
this level of relative accuracy in an hour; ‘-’ means that either the
algorithm does not achieve this level of relative accuracy in 4 hours, or it
exceeds the allocated memory.

data | $n$ | $d$ | $\log_{10}\rho$ | 2Stage-Rule2 | 2Stage-Rule4 | 2Stage-Rule1 | Rule1 | CP [6]
---|---|---|---|---|---|---|---|---
SD1 | 30000 | 4 | -3 | 6.59 (1.35) | 7.82 (2.15) | - | - | -
SD1 | 30000 | 4 | -4 | 72.78 (5.15) | 63.17 (5.42) | - | - | -
SD1 | 30000 | 10 | -3 | 39.51 (4.24) | 26.98 (2.25) | - | - | -
SD1 | 30000 | 10 | -4 | 1344.28 (158.41) | 440.54 (33.70) | - | - | -
SD1 | 30000 | 20 | -3 | 180.93 (18.49) | 96.30 (17.46) | - | - | -
SD1 | 30000 | 20 | -4 | 1831.93 (125.08) | 588.26 (48.65) | - | - | -
SD1 | 100000 | 4 | -3 | 30.83 (1.11) | 33.80 (2.55) | - | - | -
SD1 | 100000 | 4 | -4 | 105.45 (14.03) | 150.10 (20.82) | - | - | -
SD1 | 100000 | 10 | -3 | 41.94 (2.28) | 40.14 (1.89) | - | - | -
SD1 | 100000 | 10 | -4 | 846.73 (58.28) | 541.58 (81.95) | - | - | -
SD1 | 100000 | 20 | -3 | 75.02 (5.47) | 60.07 (5.91) | - | - | -
SD1 | 100000 | 20 | -4 | 2334.87 (226.08) | 1101.22 (169.24) | - | - | -
SD2 | 30000 | 4 | -3 | 8.24 (3.37) | 7.09 (0.64) | - | - | -
SD2 | 30000 | 4 | -4 | 59.41 (4.89) | 57.01 (3.68) | - | - | -
SD2 | 30000 | 10 | -3 | 45.09 (5.84) | 27.70 (1.78) | - | - | -
SD2 | 30000 | 10 | -4 | 1277.23 (114.72) | 450.98 (39.07) | - | - | -
SD2 | 30000 | 20 | -3 | 224.44 (8.14) | 116.67 (15.07) | - | - | -
SD2 | 30000 | 20 | -4 | 3018.58 (140.17) | 781.01 (58.95) | - | - | -
SD2 | 100000 | 4 | -3 | 23.48 (3.47) | 35.01 (2.82) | - | - | -
SD2 | 100000 | 4 | -4 | 106.05 (10.88) | 138.54 (18.83) | - | - | -
SD2 | 100000 | 10 | -3 | 42.74 (3.26) | 40.77 (7.13) | - | - | -
SD2 | 100000 | 10 | -4 | 820.40 (90.98) | 463.75 (34.37) | - | - | -
SD2 | 100000 | 20 | -3 | 86.63 (11.25) | 65.26 (8.63) | - | - | -
SD2 | 100000 | 20 | -4 | 2680.80 (202.38) | 1031.17 (100.00) | - | - | -
RD1 | 30000 | 4 | -4 | 18.95 (1.29) | 25.79 (1.54) | - | - | -
RD1 | 30000 | 4 | -5 | 35.16 (3.33) | 48.26 (4.82) | - | - | -
RD2 | 15000 | 4 | -4 | 51.69 (3.92) | 31.62 (2.18) | * | - | -
RD2 | 15000 | 4 | -5 | 419.94 (28.27) | 259.45 (37.73) | * | * | -

Subgradient regularization and estimation error: The $\ell_{2}$-regularization
of the subgradients can lead to better out-of-sample performance of the convex
function estimate when compared to the unregularized least squares estimate
(i.e., when $\rho=0$). This is due to the piecewise linear nature of the
convex regression fit, which can lead to overfitting near the boundary of the
convex hull of the covariates (see [2, 25]). The $\ell_{2}$-regularization on
subgradients can regulate boundary behavior and improve performance of the
estimator: see [25] for theoretical support. Here, we illustrate this
numerically for both synthetic and real datasets.

Figure 2 presents the training and test root mean squared error (RMSE), on
both real and synthetic datasets. We evaluate the RMSEs on the training and
test sets. To highlight the performance of the convex function estimate near
the boundary of the convex hull, we compute the RMSE on the boundary
points999We compute the convex hull of the training set and identify points in
the test set near the boundary of this convex hull, according to the distance
of each point to the convex hull.. Figure 2 validates our intuition that a
value of $\rho>0$ can improve the statistical performance (in terms of RMSE)
compared to the unregularized case $\rho=0$. In particular, $\rho=10^{-3}$
appears to be a good choice for the synthetic data set, and $\rho=10^{-4}\sim
10^{-5}$ for the real data set.

| RD1 with $n=1000,d=4$ | SD1 with $n=300,d=4$
---|---|---
RMSE |  |

Figure 2: Plots (in log-scale) of RMSEs evaluated on training set, test set
and boundary test set versus $\rho$’s for a real dataset (left) and synthetic
dataset (right). We consider ten replications (subsamples) and plot the median
(markers) and mean absolute deviation (error bars). We observe that the
minimum RMSE on the test/boundary set occurs when $\rho\approx
10^{-4}$-$10^{-5}$ (left panel) and $\rho=10^{-3}$ (right panel). The
unregularized fit $\rho=0$ leads to worse out-of-sample performance for both
cases.

### Acknowledgements

We would like to thank Haihao Lu for discussions in the early stages of the
work. The authors were partially supported by the Office of Naval Research
ONR-N000141812298 (Young Investigator Award), the National Science Foundation
(NSF-IIS-1718258) and MIT-IBM Watson AI Lab.

## References

  * [1] Necdet Serhat Aybat and Zi Wang.  A parallelizable dual smoothing method for large scale convex regression problems.  arXiv preprint arXiv:1608.02227, 2016.
  * [2] Fadoua Balabdaoui.  Consistent estimation of a convex density at the origin.  Mathematical Methods of Statistics, 16(2):77–95, 2007.
  * [3] Gábor Balázs.  Convex regression: Theory, practice, and applications.  PhD thesis, University of Alberta, 2016.
  * [4] Amir Beck and Marc Teboulle.  A fast iterative shrinkage-thresholding algorithm for linear inverse problems.  SIAM journal on imaging sciences, 2(1):183–202, 2009.
  * [5] Dimitri P Bertsekas.  Nonlinear programming.  Journal of the Operational Research Society, 48(3):334–334, 1997\.
  * [6] Dimitris Bertsimas and Nishanth Mundru.  Sparse convex regression.  INFORMS Journal on Computing, 33(1):262–279, 2021.
  * [7] Dimitris Bertsimas and John N Tsitsiklis.  Introduction to linear optimization, volume 6.  Athena Scientific Belmont, MA, 1997.
  * [8] Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, and Jonathan Eckstein.  Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers.  Number 3(1). Now Publishers, 2011.
  * [9] Stephen Boyd and Lieven Vandenberghe.  Convex Optimization.  Cambridge University Press, Cambridge, 2004.
  * [10] George B Dantzig and Philip Wolfe.  Decomposition principle for linear programs.  Operations research, 8(1):101–111, 1960.
  * [11] Jerome H Friedman.  Greedy function approximation: a gradient boosting machine.  Annals of statistics, pages 1189–1232, 2001.
  * [12] Avishek Ghosh, Ashwin Pananjady, Adityanand Guntuboyina, and Kannan Ramchandran.  Max-affine regression: Provable, tractable, and near-optimal statistical estimation.  arXiv preprint arXiv:1906.09255, 2019.
  * [13] P. Green and B. Silverman.  Nonparametric Regression and Generalized Linear Models: A Roughness Penalty Approach.  Chapman and Hall, London., 1994.
  * [14] Qiyang Han and Jon A Wellner.  Multivariate convex regression: global risk bounds and adaptation.  arXiv preprint arXiv:1601.06844, 2016.
  * [15] Qiyang Han and Jon A Wellner.  Convergence rates of least squares regression estimators with heavy-tailed errors.  The Annals of Statistics, 47(4):2286–2319, 2019.
  * [16] Lauren A Hannah and David B Dunson.  Multivariate convex regression with adaptive partitioning.  The Journal of Machine Learning Research, 14(1):3261–3294, 2013\.
  * [17] Alan J Hoffman.  On approximate solutions of systems of linear inequalities.  Journal of Research of the National Bureau of Standards, 49(4):263, 1952.
  * [18] Andrew L Johnson and Daniel R Jiang.  Shape constraints in economics and operations research.  Statistical Science, 33(4):527–546, 2018.
  * [19] Hamed Karimi, Julie Nutini, and Mark Schmidt.  Linear convergence of gradient and proximal-gradient methods under the polyak-łojasiewicz condition.  In Joint European Conference on Machine Learning and Knowledge Discovery in Databases, pages 795–811. Springer, 2016.
  * [20] Gil Kur, Yuval Dagan, and Alexander Rakhlin.  Optimality of maximum likelihood for log-concave density estimation and bounded convex regression.  arXiv preprint arXiv:1903.05315, 2019.
  * [21] Eunji Lim and Peter W Glynn.  Consistency of multidimensional convex regression.  Operations Research, 60(1):196–208, 2012.
  * [22] Meixia Lin, Defeng Sun, and Kim-Chuan Toh.  Efficient algorithms for multivariate shape-constrained convex regression problems.  arXiv preprint arXiv:2002.11410, 2020.
  * [23] Dong C Liu and Jorge Nocedal.  On the limited memory bfgs method for large scale optimization.  Mathematical programming, 45(1):503–528, 1989.
  * [24] Haihao Lu and Rahul Mazumder.  Randomized gradient boosting machine.  arXiv preprint arXiv:1810.10158, 2018.
  * [25] Rahul Mazumder, Arkopal Choudhury, Garud Iyengar, and Bodhisattva Sen.  A computational framework for multivariate convex regression and its variants.  Journal of the American Statistical Association, pages 1–14, 2018\.
  * [26] Maethee Mekaroonreung and Andrew L Johnson.  Estimating the shadow prices of so2 and nox for us coal power plants: a convex nonparametric least squares approach.  Energy Economics, 34(3):723–732, 2012.
  * [27] Ion Necoara, Yurii Nesterov, and Francois Glineur.  Linear convergence of first order methods for non-strongly convex optimization.  Mathematical Programming, pages 1–39, 2018.
  * [28] Yu Nesterov.  Gradient methods for minimizing composite functions.  Mathematical Programming, 140(1):125–161, 2013.
  * [29] Yurii Nesterov.  Lectures on convex optimization.  Springer, 2018.
  * [30] Jorge Nocedal and Stephen Wright.  Numerical optimization.  Springer Science & Business Media, 2006.
  * [31] Brendan O’donoghue and Emmanuel Candes.  Adaptive restart for accelerated gradient schemes.  Foundations of computational mathematics, 15(3):715–732, 2015.
  * [32] Fred Ramsey and Daniel Schafer.  The statistical sleuth: a course in methods of data analysis.  Cengage Learning, 2012.
  * [33] T. Robertson, F.T. Wright, and R. Dykstra.  Order Restricted Statistical Inference.  Probability and Statistics Series. Wiley, 1988.
  * [34] Richard Samworth and Bodhisattva Sen.  Editorial: Special issue on “nonparametric inference under shape constraints”.  Statistical Science, 2018.
  * [35] Emilio Seijo and Bodhisattva Sen.  Nonparametric least squares estimation of a multivariate convex regression function.  The Annals of Statistics, 39(3):1633–1657, 2011.
  * [36] Alexander Shapiro, Darinka Dentcheva, and Andrzej Ruszczyński.  Lectures on stochastic programming: modeling and theory.  SIAM, 2009.
  * [37] Hal R Varian.  The nonparametric approach to demand analysis.  Econometrica: Journal of the Econometric Society, pages 945–973, 1982.
  * [38] Grace Wahba.  Spline models for observational data, volume 59.  Siam, 1990.
  * [39] Yongqiao Wang and Shouyang Wang.  Estimating $\alpha$-frontier technical efficiency with shape-restricted kernel quantile regression.  Neurocomputing, 101:243–251, 2013.

## Appendix and Supplementary Material

##  Appendix A Proofs

In this section, we present the proofs of Lemma 1 and Theorem 1. The proof of
Theorem 1 is based on Lemmas 1 and 4. We first present the proof of Lemma 1
(Sec. A.1); then present Lemma 4 (Sec. A.2), followed by the proof of Theorem
1.

###  A.1 Proof of Lemma 1

For any $\bm{\theta}\in\mathbb{R}^{N}$ and $P\leq N$, we define two auxiliary
norms:

$\lVert\bm{\theta}\rVert_{\mathcal{G}[P]}=\left(\sum_{i\in\mathcal{G}[[N],P,\bm{\theta}]}\theta_{i}^{2}\right)^{1/2}~{}~{}~{}\text{and}~{}~{}~{}\lVert\bm{\theta}\rVert_{\mathcal{U}[P]}=\mathbb{E}\left(\sum_{i\in\mathcal{U}[[N],P]}\theta_{i}^{2}\right)^{1/2}$

where, $\mathcal{G}[[N],P,\bm{\theta}]$, $\mathcal{U}[P]$ are defined in
Definition 1 and $\mathbb{E}(\cdot)$ is the expectation wrt scheme
${\mathcal{U}}$. Proposition 3 links these norms to the Euclidean norm:

######  Proposition 3.

The following statement holds:

$\lVert\bm{\theta}\rVert^{2}\geq\lVert\bm{\theta}\rVert_{\mathcal{G}[P]}^{2}\geq\frac{P}{N}\lVert\bm{\theta}\rVert^{2}=\lVert\bm{\theta}\rVert_{\mathcal{U}[P]}^{2}.$

###### Proof of Proposition 3.

It is easy to see that
$\lVert\bm{\theta}\rVert_{\mathcal{U}[P]}^{2}=(P/N)\lVert\bm{\theta}\rVert^{2}.$
Notice that

$\displaystyle\lVert\bm{\theta}\rVert_{\mathcal{G}[P]}^{2}~{}=~{}\max_{\bm{\pi}}\sum_{\omega\in[N]}\pi_{\omega}|\theta_{\omega}|^{2}~{}~{}~{}\text{s.t.}~{}~{}~{}\sum_{\omega\in[N]}\pi_{\omega}\leq
P,~{}~{}~{}0\leq\pi_{\omega}\leq 1,\forall\omega.$

Since $\pi_{\omega}={P}/{N}~{}\forall\omega$, is feasible for the above
problem, it follows that

$\lVert\bm{\theta}\rVert_{\mathcal{G}[P]}^{2}\geq\frac{P}{N}\lVert\bm{\theta}\rVert^{2}=\lVert\bm{\theta}\rVert_{\mathcal{U}[P]}^{2}.$

Equality above is attained if and only if $\theta_{\omega}=C$ $\forall\omega$
for some $C$. Furthermore, we note that

$\lVert\bm{\theta}\rVert_{\mathcal{G}[P]}^{2}=\sum_{\omega\in\mathcal{G}[[N],P]}|\theta_{\omega}|^{2}\leq\sum_{\omega\in[N]}|\theta_{\omega}|^{2}=\lVert\bm{\theta}\rVert^{2},$

and this equality is attained if and only if $\theta_{\omega}=0$ for all
$\omega\not\in\mathcal{G}[[N],P]$, i.e. the $N-P$ smallest values of
$|\theta_{\omega}|$ are 0. ∎

###### Proof of Lemma 1.

We divide the proof into 5 parts depending upon the 5 rules.

Rule 1: Greedy within each Block.
($\delta_{1}(\bm{\theta},\\{\Omega_{i}\\})=\cup_{i=1}^{n}\mathcal{G}[\Omega_{i},P,\bm{\theta}]$,
$\alpha_{\\{1\\}}=\frac{n-1}{P}$, and $\beta_{\\{1\\}}=1$.) For this selection
rule, we have

$\displaystyle\lVert\bm{\theta}\rVert_{\\{1\\}}^{2}$
$\displaystyle=\sum_{i=1}^{n}\lVert\bm{\theta}_{\Omega_{i}}\rVert_{\mathcal{G}[P]}^{2}.$

It is easy to see that $\|{\bm{\theta}}\|_{\\{1\\}}$ is a norm. It follows
from Proposition 3 that

$\lVert\bm{\theta}\rVert^{2}=\sum_{i=1}^{n}\lVert\bm{\theta}_{\Omega_{i}}\rVert^{2}\geq\sum_{i=1}^{n}\lVert\bm{\theta}_{\Omega_{i}}\rVert_{\mathcal{G}[P]}^{2}\geq\sum_{i=1}^{n}\frac{n-1}{P}\lVert\bm{\theta}_{\Omega_{i}}\rVert^{2}=\frac{n-1}{P}\lVert\bm{\theta}\rVert^{2}.$

Therefore, we have $\alpha_{\\{1\\}}=\frac{n-1}{P}$, and $\beta_{\\{1\\}}=1.$

Rule 2: Random.
($\delta_{2}(\bm{\theta},\\{\Omega_{i}\\})=\mathcal{U}[\Omega,K]$,
$\alpha_{\\{2\\}}=\frac{n(n-1)}{K}$ and $\beta_{\\{2\\}}=\frac{n(n-1)}{K}$).

For this selection rule, we have

$\lVert\bm{\theta}\rVert_{\\{2\\}}^{2}=\lVert\bm{\theta}\rVert_{\mathcal{U}[K]}^{2}=\frac{n(n-1)}{K}\lVert\bm{\theta}\rVert^{2}.$

Thus, in this case the norm-equivalence constants are
$\alpha_{\\{2\\}}=\beta_{\\{2\\}}=\frac{n(n-1)}{K}$.

Rule 3: Random within each Block.
($\delta_{3}(\bm{\theta},\\{\Omega_{i}\\})=\cup_{i=1}^{n}\mathcal{U}[\Omega_{i},P]$,
$\alpha_{\\{3\\}}=\frac{n-1}{P}$, and $\beta_{\\{3\\}}=\frac{n-1}{P}$.) For
this selection rule, we have

$\lVert\bm{\theta}\rVert_{\\{3\\}}^{2}=\sum_{i=1}^{n}\lVert\bm{\theta}_{\Omega_{i}}\rVert_{\mathcal{U}[P]}^{2}=\sum_{i=1}^{n}\frac{n-1}{P}\lVert\bm{\theta}_{\Omega_{i}}\rVert^{2}=\frac{n-1}{P}\lVert\bm{\theta}\rVert^{2}.$

Hence, it follows that $\alpha_{\\{3\\}}=\beta_{\\{3\\}}=\frac{n-1}{P}.$

Rule 4: Random then Greedy.
($\delta_{4}(\bm{\theta},\\{\Omega_{i}\\})=\mathcal{G}[\mathcal{U}[\Omega,M],K,\bm{\theta}]$,
$\alpha_{\\{4\\}}=\frac{n(n-1)}{K}$, and $\beta_{\\{4\\}}=\frac{n(n-1)}{M}$.)
We adapt the proof of [24] to show that under this selection rule

$\lVert\bm{\theta}\rVert_{\\{4\\}}^{2}=\sum_{l=1}^{n(n-1)}\pi(l)|\theta_{(l)}|^{2},$

where $|\theta_{(l)}|$ is the $l$-th largest value in
$\\{|\theta_{\omega}|\\}_{\omega\in\Omega}$ and $\pi(l)$ is given by

$\pi(l)=\frac{M}{n(n-1)}\sum_{k=1}^{K}\frac{\binom{l-1}{k-1}\binom{n(n-1)-l}{M-k}}{\binom{n(n-1)-1}{M-1}}.$

Here, by convention, we define $\binom{N}{\alpha}=0$ if $\alpha<0$ or
$\alpha>N$.

Let $\pi(l)$ be the probability that $|\theta_{(l)}|$ is selected. Since the
subsample is selected uniformly at random, it suffices to count the number of
combinations that include $|\theta_{(l)}|$ and those in which $|\theta_{(l)}|$
ranks among the top $K$ values. This is equivalent to choosing $k-1(\leq K-1)$
elements from $\\{|\theta_{(s)}|\\}_{s\leq l-1}$, selecting the element
$|\theta_{(l)}|$ and then choosing the remaining $(M-k)$ elements from the
rest. Therefore, the number of such combinations is

$N(l)=\sum_{k=1}^{K}\binom{l-1}{k-1}\binom{n(n-1)-l}{M-k}.$

Thus,

$\pi(l)=\frac{N(l)}{\binom{n(n-1)}{M}}=\frac{M}{n(n-1)}\sum_{k=1}^{K}\frac{\binom{l-1}{k-1}\binom{n(n-1)-l}{M-k}}{\binom{n(n-1)-1}{M-1}}.$

Notice that when $l\leq K$ (i.e., $l-1\leq K-1$), then

$N(l)=\sum_{k=1}^{l}\binom{l-1}{k-1}\binom{n(n-1)-l}{M-k}=\binom{n(n-1)-1}{M-1},$

and thus $\pi(l)=\frac{M}{n(n-1)}.$

When $l\geq n(n-1)-(M-K)+1$ and $M-k\geq M-1-K>n(n-1)-l$, then $N(l)=0$ and
thus $\pi(l)=0$.

As each element appears in the same number of combinations of size $M$ (in the
random selection step), and the greedy selection step favors the larger one,
we have the following ordering: $\pi(1)\geq\ldots\geq\pi(n(n-1))$. Therefore,
$\lVert\bm{\theta}\rVert_{\\{4\\}}$ is a norm.

Since
$\pi(l)=\mathbb{E}[\mathbb{I}_{\\{(l)\in\mathcal{G}[\mathcal{U}[\Omega,M],K]\\}}]$,
it follows that

$\sum_{l=1}^{n(n-1)}\pi(l)=\mathbb{E}\left[\sum_{l=1}^{n(n-1)}\mathbb{I}_{\\{(l)\in\mathcal{G}[\mathcal{U}[\Omega,M],K]\\}}\right]=K.$

Then, we have

$\lVert\bm{\theta}\rVert_{\\{4\\}}^{2}=\sum_{l=1}^{n(n-1)}\pi(l)|\theta_{(l)}|^{2}\geq\sum_{l=1}^{n(n-1)}\frac{K}{n(n-1)}|\theta_{(l)}|^{2}=\frac{K}{n(n-1)}\lVert\bm{\theta}\rVert^{2}.$
(10)

On the other hand, since $\pi(l)\leq\pi(1)=\frac{M}{n(n-1)}$, it follows that

$\lVert\bm{\theta}\rVert_{\\{4\\}}^{2}=\sum_{l=1}^{n(n-1)}\pi(l)|\theta_{(l)}|^{2}\leq\sum_{l=1}^{n(n-1)}\frac{M}{n(n-1)}|\theta_{(l)}|^{2}=\frac{M}{n(n-1)}\lVert\bm{\theta}\rVert^{2}.$
(11)

Hence from (10) and (11), we obtain $\alpha_{\\{4\\}}=\frac{n(n-1)}{K}$ and
$\beta_{\\{4\\}}=\frac{n(n-1)}{M}$.

Rule 5: Random Blocks then Greedy within each Block.
($\delta_{5}(\bm{\theta},\\{\Omega_{i}\\})=\cup_{i\in\mathcal{U}[[n],G]}\mathcal{G}[\Omega_{i},P,\bm{\theta}]$,
$\alpha_{\\{5\\}}=\frac{n(n-1)}{GP}$, and $\beta_{\\{5\\}}=\frac{n}{G}$.)

By the selection rule, we have

$\lVert\bm{\theta}\rVert_{\\{5\\}}^{2}=\frac{G}{n}\sum_{i=1}^{n}\lVert\bm{\theta}_{\Omega_{i}}\rVert_{\mathcal{G}[P]}^{2}=\frac{G}{n}\lVert\bm{\theta}\rVert_{\\{1\\}}^{2}$

with $\alpha_{\\{5\\}}=\frac{n}{G}\alpha_{\\{1\\}}=\frac{n(n-1)}{GP}$ and
$\beta_{\\{5\\}}=\frac{n}{G}\beta_{\\{1\\}}=\frac{n}{G}.$

∎

###  A.2 Auxiliary lemmas for the proof of Theorem 1

Here we present Lemmas 3 and 4 that will be used for the proof of Theorem 1.

######  Lemma 2 (Hoffman’s Lemma [17]).

Let $\Lambda=\\{\bm{\lambda}:\bm{D\lambda}=\bm{s},\bm{\lambda}\leq\bm{0}\\}$.
There is a constant $\mu_{\bm{D}}>0$ that depends only on $\bm{D}$ such that
for any $\bm{\lambda}\leq\bm{0}$, there exists an $\bm{\lambda}_{0}\in\Lambda$
with

$\lVert\bm{D\lambda}-\bm{s}\rVert^{2}\geq\mu_{\bm{D}}\lVert\bm{\lambda}-\bm{\lambda}_{0}\rVert^{2}.$

The constant $\mu_{\bm{D}}$ is called the Hoffman’s constant associated with
$\bm{D}$.

Let $\bm{D}^{\star}=\begin{bmatrix}\bm{A}^{\top}\\\
\bm{B}^{\top}/\sqrt{\rho}\end{bmatrix}$ be the coefficient matrix
corresponding to set of optimal solutions $\bm{\Lambda}^{\star}$ with
$W=\Omega$ in (3). We get the following sufficient descent lemma according to
[19]:

######  Lemma 3 ([19] Case 3, Page 18).

For any $\bm{\lambda}\leq\bm{0}$ and optimal dual solution
$\bm{\lambda}^{\star}$ of ($D$), we have

$L(\bm{\lambda}^{\star})\geq
L(\bm{\lambda})-\frac{1}{2\mu_{\bm{D}^{\star}}}\mathcal{D}(\bm{\lambda},\nabla
L(\bm{\lambda}),\mu_{\bm{D}^{\star}}),$ (12)

where $\mu_{\bm{D}^{\star}}>0$ is the Hoffman’s constant associated with
$\bm{D}^{\star}$.

######  Definition 4.

Define a function
$\mathcal{D}(\cdot,\cdot,\cdot):\mathbb{R}_{-}^{N}\times\mathbb{R}^{N}\times\mathbb{R}_{++}\to\mathbb{R}_{+}$
as

$\mathcal{D}(\bm{\lambda},\bm{\theta},\mu)=-2\mu\min_{\bm{\lambda}^{\prime}\leq
0}\left\\{\langle\bm{\theta},\bm{\lambda}^{\prime}-\bm{\lambda}\rangle+\frac{\mu}{2}\lVert\bm{\lambda}^{\prime}-\bm{\lambda}\rVert^{2}\right\\}.$

Lemma 4 makes use of the following useful proposition.

######  Proposition 4.

Define
$d(\cdot,\cdot,\cdot):\mathbb{R}_{-}\times\mathbb{R}_{+}\times\mathbb{R}_{++}\to\mathbb{R}$
as

$d(\lambda,\theta,\mu)=-2\mu\min_{\lambda^{\prime}\leq
0}\left\\{\theta(\lambda^{\prime}-\lambda)+\frac{\mu}{2}(\lambda^{\prime}-\lambda)^{2}\right\\}.$

We have the following:

  * (a)

if $\lambda=0$, then for any $\mu>0$,
$d(\lambda,\theta,\mu)=\max\\{\theta,0\\}^{2};$

  * (b)

if $\lambda\theta=0$ with $\lambda\leq 0$ and $\theta\leq 0$, then for any
$\mu>0$, $d(\lambda,\theta,\mu)=0$;

  * (c)

if $\lambda\leq 0$, then for any $\mu>0$, $d(\lambda,\theta,\mu)\geq 0$.

######  Proof..

Part (a): If $\lambda=0$, then

$d(\lambda,\theta,\mu)=-2\mu\min_{\lambda^{\prime}\leq
0}\left\\{\theta\lambda^{\prime}+\frac{\mu}{2}(\lambda^{\prime})^{2}\right\\}=\left\\{\begin{array}[]{ll}0&\text{
if }\theta\leq 0\\\ \theta^{2}&\text{ if }\theta>0\end{array}\right.,$

i.e. $d(\lambda,\theta,\mu)=\max\\{\theta,0\\}^{2}$.

Part (b): If $\lambda\theta=0$ with $\lambda=0$ and $\theta\leq 0$, then
$d(\lambda,\theta,\mu)=0$ follows from Part (a).

If $\lambda\theta=0$ with $\lambda<0$, then $\theta=0$, so
$d(\lambda,\theta,\mu)=-2\mu\min_{\lambda^{\prime}\leq
0}\\{\frac{\mu}{2}(\lambda^{\prime}-\lambda)^{2}\\}=0.$

Part (c): If $\lambda\leq 0$, then $\lambda^{\prime}=\lambda\leq 0$ is always
feasible, so $d(\lambda,\theta,\mu)\geq 0.$ ∎

######  Lemma 4.

Suppose there exists $\alpha\geq 1$ such that for any $m$ and $W^{0}$ the
following

$\alpha\mathbb{E}_{\eta_{m+1}}\left[\sum_{\omega\in W^{m+1}\backslash
W^{m}}\max\\{\nabla_{\lambda_{\omega}}L(\bm{\lambda}^{m}),0\\}^{2}\right]\geq\sum_{\omega\in\Omega\backslash
W^{m}}\max\\{\nabla_{\lambda_{\omega}}L(\bm{\lambda}^{m}),0\\}^{2}$ (13)

holds almost surely, where $\mathbb{E}_{\eta_{m+1}}$ denotes expectation over
all the random sources in $(m+1)$-th iteration (conditional on events up to
iteration $m$). Then, the following holds:

$\mathbb{E}_{\eta_{m+1}}[L(\bm{\lambda}^{m+1})-L^{\star}]\leq\left(1-\frac{\mu_{\bm{D}^{\star}}}{\alpha\sigma}\right)[L(\bm{\lambda}^{m})-L^{\star}].$

###### Proof of Lemma 4.

We make use of the following notation in the proof:

$\text{$\lambda_{\omega}^{m}$ is the $\omega$-th component of
$\bm{\lambda}^{m}$},~{}~{}~{}~{}~{}\nabla^{m}_{\omega}=\nabla_{\lambda_{\omega}}L(\bm{\lambda}^{m})~{}~{}~{}\text{and}~{}~{}~{}~{}\tilde{\nabla}^{m}_{\omega}=\max\\{\nabla^{m}_{\omega},0\\}.$

The proof has three steps, and we consider both exact and inexact cases in
each step.

Step 1: In the first part, we will show that

$\displaystyle
L(\bm{\lambda}^{m})-L(\bm{\lambda}^{\star})~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\leq$
$\displaystyle\left\\{\begin{array}[]{ll}\displaystyle\frac{1}{2\mu_{\bm{D}^{\star}}}\sum_{\omega\in\Omega\backslash
W^{m}}(\tilde{\nabla}_{\omega}^{m})^{2}&\text{ (exact case)}\\\
\displaystyle\frac{1}{2\mu_{\bm{D}^{\star}}}\left[\sum_{\omega\in
W^{m}}d(\lambda^{m}_{\omega},\nabla^{m}_{\omega},\mu_{\bm{D}^{\star}})+\sum_{\omega\in\Omega\backslash
W^{m}}(\tilde{\nabla}_{\omega}^{m})^{2}\right]&\text{ (inexact
case)}.\end{array}\right.$ (16)

It follows from Lemma 3 that

$L(\bm{\lambda}^{\star})\geq
L(\bm{\lambda}^{m})-\frac{1}{2\mu_{\bm{D}^{\star}}}\mathcal{D}(\bm{\lambda}^{m},\nabla
L(\bm{\lambda}^{m}),\mu_{\bm{D}^{\star}}).$ (17)

By the definitions of ${\mathcal{D}}(\cdot)$ (Definition 4) and $d(\cdot)$
(Proposition 4), we have

$\displaystyle\mathcal{D}(\bm{\lambda}^{m},\nabla
L(\bm{\lambda}^{m}),\mu_{\bm{D}^{\star}})$
$\displaystyle=\sum_{\omega\in\Omega}d(\lambda^{m}_{\omega},\nabla^{m}_{\omega},\mu_{\bm{D}^{\star}})$
$\displaystyle=\sum_{\omega\in
W^{m}}d(\lambda^{m}_{\omega},\nabla^{m}_{\omega},\mu_{\bm{D}^{\star}})+\sum_{\omega\in\Omega\backslash
W^{m}}d(\lambda^{m}_{\omega},\nabla^{m}_{\omega},\mu_{\bm{D}^{\star}})$ (18)

##### Inexact Case:

By the definition of $\bm{\lambda}^{W^{m}}$, we know that
$\lambda_{\omega}^{m}=0$ for $\omega\not\in W^{m}$. Then it follows from
Proposition 4 that

$\displaystyle\mathcal{D}(\bm{\lambda}^{m},\nabla
L(\bm{\lambda}^{m}),\mu_{\bm{D}^{\star}})=\sum_{\omega\in
W^{m}}d(\lambda^{m}_{\omega},\nabla^{m}_{\omega},\mu_{\bm{D}^{\star}})+\sum_{\omega\in\Omega\backslash
W^{m}}(\tilde{\nabla}_{\omega}^{m})^{2}$ (19)

##### Exact Case:

By the primal feasibility, dual feasibility and complementary slackness of
$\bm{\lambda}_{W^{m}}^{\star}$, we have
$\lambda_{\omega}^{m}\nabla_{\omega}^{m}=0$ with
$\lambda_{\omega}^{m},\nabla_{\omega}^{m}\leq 0$ for $\omega\in W^{m}$.
Combining this with $\lambda_{\omega}^{m}=0$ for $\omega\not\in W^{m}$, by
Proposition 4, we have

$\mathcal{D}(\bm{\lambda}^{m},\nabla
L(\bm{\lambda}^{m}),\mu_{\bm{D}^{\star}})=\sum_{\omega\in\Omega\backslash
W^{m}}(\tilde{\nabla}_{\omega}^{m})^{2}.$ (20)

The conclusions from (19) and (20) lead to (A.2).

Step 2: In the second step, we will show that

$\displaystyle
L(\bm{\lambda}^{m+1})-L(\bm{\lambda}^{m})~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\leq$
$\displaystyle\left\\{\begin{array}[]{ll}\displaystyle-\frac{1}{2\sigma}\sum_{\omega\in
W^{m+1}\backslash W^{m}}(\tilde{\nabla}_{\omega}^{m})^{2}&\text{ (exact
case)}\\\ \displaystyle-\frac{1}{2\sigma}\left[\sum_{\omega\in
W^{m}}d(\lambda^{m}_{\omega},\nabla^{m}_{\omega},\sigma)+\sum_{\omega\in
W^{m+1}\backslash W^{m}}(\tilde{\nabla}_{\omega}^{m})^{2}\right]&\text{
(inexact case)}\end{array}\right.$ (23)

##### Exact Case:

Let
$\Lambda^{m+1}=\\{\bm{\lambda}\in\mathbb{R}^{n(n-1)}:\lambda_{\omega}=0,\forall\omega\in\Omega\backslash
W^{m+1}\\}.$ Recall that $\bm{\lambda}^{m+1}$ minimizes $L(\bm{\lambda})$ over
$\Lambda^{m+1}$. Therefore, we have the following:

$\displaystyle L(\bm{\lambda}^{m+1})$
$\displaystyle=\min_{\bm{\lambda}\in\Lambda^{m+1}}L(\bm{\lambda})$
$\displaystyle\leq\min_{\bm{\lambda}\in\Lambda^{m+1}}\left\\{L(\bm{\lambda}^{m})+\langle\nabla
L(\bm{\lambda}^{m}),\bm{\lambda}-\bm{\lambda}^{m}\rangle+\frac{\sigma}{2}\lVert\bm{\lambda}-\bm{\lambda}^{m}\rVert^{2}\right\\}$
$\displaystyle=L(\bm{\lambda}^{m})-\frac{1}{2\sigma}\sum_{\omega\in
W^{m}}d(\lambda^{m}_{\omega},\nabla^{m}_{\omega},\sigma)-\frac{1}{2\sigma}\sum_{\omega\in
W^{m+1}\backslash W^{m}}d(\lambda^{m}_{\omega},\nabla^{m}_{\omega},\sigma)$
$\displaystyle=L(\bm{\lambda}^{m})-\frac{1}{2\sigma}\sum_{\omega\in
W^{m+1}\backslash W^{m}}(\tilde{\nabla}_{\omega}^{m})^{2}$ (24)

where the first inequality uses $\sigma$-smoothness of $L$; the last line
follows from an argument similar to (20) where we use Proposition 4 and
complementary slackness.

##### Inexact Case:

Here we take one (or more) projected gradient step(s) to partially minimize
the reduced dual. Let $\bm{\lambda}^{m+\frac{1}{2}}$ be obtained after one
projected gradient descent step from $\bm{\lambda}^{m}$. Then we have

$\displaystyle L(\bm{\lambda}^{m+1})$ $\displaystyle\leq
L(\bm{\lambda}^{m+\frac{1}{2}})$ $\displaystyle\leq
L(\bm{\lambda}^{m})+\langle\nabla
L(\bm{\lambda}^{m}),\bm{\lambda}^{m+\frac{1}{2}}-\bm{\lambda}^{m}\rangle+\frac{\sigma}{2}\lVert\bm{\lambda}^{m+\frac{1}{2}}-\bm{\lambda}^{m}\rVert^{2}$
$\displaystyle=L(\bm{\lambda}^{m})-\frac{1}{2\sigma}\sum_{\omega\in
W^{m+1}}d(\lambda^{m}_{\omega},\nabla^{m}_{\omega},\sigma)$
$\displaystyle=L(\bm{\lambda}^{m})-\frac{1}{2\sigma}\sum_{\omega\in
W^{m}}d(\lambda^{m}_{\omega},\nabla^{m}_{\omega},\sigma)-\frac{1}{2\sigma}\sum_{\omega\in
W^{m+1}\backslash W^{m}}(\tilde{\nabla}_{\omega}^{m})^{2},$ (25)

where the first inequality follows from the descent property of projected
gradient descent; the second inequality uses $\sigma$-smoothness of $L$; and
the last line follows from Proposition 4.

Finally, the result in (A.2) follows from (24) and (25)

Step 3: For the third step, we will show that the following holds

$\sigma\alpha\mathbb{E}_{\eta_{m+1}}[L(\bm{\lambda}^{m+1})-L(\bm{\lambda}^{m})]\leq-\mu_{\bm{D}^{\star}}[L(\bm{\lambda}^{m})-L(\bm{\lambda}^{\star})].$
(26)

##### Exact Case:

For this case, we have the following chain of inequalities:

$\displaystyle-2\sigma\alpha\mathbb{E}_{\eta_{m+1}}[L(\bm{\lambda}^{m+1})-L(\bm{\lambda}^{m})]$
$\displaystyle\geq\alpha\mathbb{E}_{\eta_{m+1}}\sum_{\omega\in
W^{m+1}\backslash W^{m}}(\tilde{\nabla}_{\omega}^{m})^{2}$
$\displaystyle\geq\sum_{\omega\in\Omega\backslash
W^{m}}(\tilde{\nabla}_{\omega}^{m})^{2}$ $\displaystyle\geq
2\mu_{\bm{D}^{\star}}[L(\bm{\lambda}^{m})-L(\bm{\lambda}^{\star})],$ (27)

where the first inequality uses (A.2); the second inequality is the assumption
(13); and the last line uses (A.2).

##### Inexact Case:

Since $\lambda_{\omega}^{m}\leq 0$, we know
$d(\lambda_{\omega}^{m},\nabla_{\omega}^{m},\sigma)\geq 0$ by Proposition 4.
Using the fact that $\alpha\geq 1$, we have the following:

$\displaystyle-2\sigma\alpha\mathbb{E}_{\eta_{m+1}}[L(\bm{\lambda}^{m+1})-L(\bm{\lambda}^{m})]$
$\displaystyle\geq\alpha\sum_{\omega\in
W^{m}}d(\lambda^{m}_{\omega},\nabla^{m}_{\omega},\sigma)+\alpha\mathbb{E}_{\eta_{m+1}}\sum_{\omega\in
W^{m+1}\backslash W^{m}}(\tilde{\nabla}_{\omega}^{m})^{2}$
$\displaystyle\geq\sum_{\omega\in
W^{m}}d(\lambda^{m}_{\omega},\nabla^{m}_{\omega},\sigma)+\sum_{\omega\in\Omega\backslash
W^{m}}(\tilde{\nabla}_{\omega}^{m})^{2}$ $\displaystyle\geq
2\mu_{\bm{D}^{\star}}[L(\bm{\lambda}^{m})-L(\bm{\lambda}^{\star})],$ (28)

where the first inequality uses (A.2); the second inequality uses assumption
(13), $\alpha\geq 1$ and
$d(\lambda_{\omega}^{m},\nabla_{\omega}^{m},\sigma)\geq 0$; and the last line
uses (A.2).

The statement in (26) follows from (27) and (28).

Finally, we complete the proof by using (26) and observing that:

$\displaystyle\mathbb{E}_{\eta_{m+1}}[L(\bm{\lambda}^{m+1})-L^{\star}]=$
$\displaystyle\mathbb{E}_{\eta_{m+1}}[L(\bm{\lambda}^{m+1})-L(\bm{\lambda}^{m})]+L(\bm{\lambda}^{m})-L^{\star}$
$\displaystyle\leq$
$\displaystyle\left(1-\frac{\mu_{\bm{D}^{\star}}}{\alpha\sigma}\right)[L(\bm{\lambda}^{m})-L^{\star}].$

∎

###  A.3 Proof of Theorem 1

The proof of Theorem 1 uses Lemmas 1 and 4.

###### Proof of Theorem 1.

Recall that $\lambda_{\omega}^{m}$ is the $\omega$-th component of
$\bm{\lambda}^{m}$,
$\nabla^{m}_{\omega}=\nabla_{\lambda_{\omega}}L(\bm{\lambda}^{m})$, and
$\tilde{\nabla}^{m}_{\omega}=\max\\{\nabla^{m}_{\omega},0\\}$. Now let
$\\{\Omega_{i}\\}$ be one of the partitions $\\{\Omega_{i\cdot}\\}$ or
$\\{\Omega_{\cdot i}\\}$, and $W_{i}^{m}$ be the corresponding partition for
$W^{m}$. Let $\Delta=W^{m+1}\backslash W^{m}$. Using this notation, the
condition of Lemma 4 reduces to

$\alpha\mathbb{E}_{\eta_{m+1}}[\sum_{\omega\in\Delta}(\tilde{\nabla}_{\omega}^{m})^{2}]\geq\sum_{\omega\in\Omega\backslash
W^{m}}(\tilde{\nabla}_{\omega}^{m})^{2}.$ (29)

We organize the proof into four steps: (1) we construct a random vector
$\bm{g}\in\mathbb{R}^{n(n-1)}$ from $\bm{\nabla}^{m}=\nabla
L(\bm{\lambda}^{m})$; (2) we then show $\sum_{\omega\in\Omega\backslash
W^{m}}(\tilde{\nabla}_{\omega}^{m})^{2}$ equals $\lVert\bm{g}\rVert^{2}$; (3)
we relate
$\mathbb{E}_{\eta_{m+1}}[\sum_{\omega\in\Delta}(\tilde{\nabla}_{\omega}^{m})^{2}]$
to $\lVert\bm{g}\rVert_{\\{\ell\\}}^{2}$; and (4) finally, we apply Lemmas 1
and 4 to complete the proof.

Step 1: Define each entry $g_{\omega}$ of $\bm{g}$ as follows:

$g_{\omega}=\left\\{\begin{array}[]{ll}\tilde{\nabla}_{\omega}^{m}&\text{ if
}\omega\in\Omega\backslash W^{m}\\\ 0&\text{ if }\omega\in
W^{m}\end{array}\right..$ (30)

Notice $\bm{g}$ is a random vector depending upon random sources from
iterations $1,\ldots,m$.

Step 2: By the definition of $\bm{g}$, we have
$\sum_{\omega\in\Omega\backslash
W^{m}}(\tilde{\nabla}_{\omega}^{m})^{2}=\lVert\bm{g}\rVert_{2}^{2}.$

Step 3: Recall that $\Delta^{\prime}$ denotes the set of selected pairs as per
a rule, and we consider a subset
$\Delta=\\{\omega\in\Delta^{\prime}:\tilde{\nabla}_{\omega}^{m}=\nabla_{\omega}^{m}>0\\}$—that
is, $\tilde{\nabla}_{\omega}^{m}=0$ for any
$\omega\in\Delta^{\prime}\backslash\Delta$. Thus,

$\mathbb{E}_{\eta_{m+1}}\left[\sum_{\omega\in\Delta}(\tilde{\nabla}_{\omega}^{m})^{2}\right]=\mathbb{E}_{\eta_{m+1}}\left[\sum_{\omega\in\Delta^{\prime}}(\tilde{\nabla}_{\omega}^{m})^{2}\right].$

Note that
$\Delta^{\prime}=\delta_{\ell}(\tilde{\bm{\nabla}}^{m}_{\Omega\backslash
W^{m}},\\{\Omega_{i}\setminus W_{i}^{m}\\})$, and let
$\Delta^{\prime\prime}=\delta_{\ell}(\bm{g},\\{\Omega_{i}\\})$. Notice that
$\bm{g}$ has more zero coordinates compared to $\tilde{\bm{\nabla}}^{m}$ (see
(30)). Thus, we have

$\mathbb{E}_{\eta_{m+1}}\left[\sum_{\omega\in\Delta}(\tilde{\nabla}_{\omega}^{m})^{2}\right]=\mathbb{E}_{\eta_{m+1}}\left[\sum_{\omega\in\Delta^{\prime}}(\tilde{\nabla}_{\omega}^{m})^{2}\right]\geq\mathbb{E}_{\eta_{m+1}}\left[\sum_{\omega\in\Delta^{\prime\prime}}(g_{\omega})^{2}\right]=\lVert\bm{g}\rVert_{\\{\ell\\}}^{2},$

where the last equality follows from Definition 3.

Step 4: From Step 2, Step 3 and Lemma 1, we arrive at (29). Therefore, it
follows from Lemma 4 that

$\mathbb{E}_{\eta_{m+1}}[L(\bm{\lambda}^{m})-L^{\star}]\leq\left(1-\frac{\mu_{\bm{D}^{\star}}}{\alpha_{\\{\ell\\}}\sigma}\right)[L(\bm{\lambda}^{m})-L^{\star}]$

and (using tower property of expectation) we arrive at the conclusion of the
theorem:

$\mathbb{E}[L(\bm{\lambda}^{m})-L^{\star}]\leq\left(1-\frac{\mu_{\bm{D}^{\star}}}{\alpha_{\\{\ell\\}}\sigma}\right)^{m}[L(\bm{\lambda}^{0})-L^{\star}].$

∎

######  Remark 4.

Both accelerated gradient methods (APG) and L-BFGS can be used to solve the
reduced problems to optimality—so, the proof for the exact subproblem
optimization (in Theorem 1) for these cases will be the same as that for PGD,
The proof for the inexact subproblem optimization for PGD uses only the
sufficient decrease condition of the first PGD update and the fact that PGD is
a descent algorithm. Since APG with adaptive restart (function scheme) [31]
and L-BFGS [23] are descent algorithms, the theory presented above, goes
through as long as the progress made by APG or L-BFGS is at least as large as
the progress made by the first step of APG. To this end, for L-BFGS, when
setting the initial inverse Hessian approximation as the identity matrix, i.e.
$\bm{H}_{0}=\bm{I}$, the first step of L-BFGS is exactly PGD update, so the
theory works for L-BFGS. For APG with restarts, if we perform a PGD update
preceding the APG updates, then the theory will go through as well.

##  Appendix B Additional Technical Details

###  B.1 Examples of unavoidable factor $\alpha_{\\{\ell\\}}$

For simplicity, we consider an unconstrained problem with objective
$f(x)=x^{\top}Hx$, where
$H=\mathrm{diag}(\sigma,1,1,\ldots,1)\in\mathbb{R}^{p\times p}$ with
$\sigma>1$. Then, $f$ is $1$-strongly convex and $\sigma$-smooth, and thus the
step size will be $1/\sigma$.

Gradient Descent: Starting from any $x^{k}=(x_{i}^{k}:i\in[p])$, the gradient
step yields

$x^{k+1}=(0,\gamma x^{k}_{2},\ldots,\gamma x^{k}_{m}),$

where $\gamma=1-\frac{1}{\sigma}$. Therefore, for $k\geq 1$,
$x^{k+1}_{1}=x^{k}_{1}=0$, and $f(x^{k+1})=\gamma^{2}f(x^{k})$. Since the
minimum objective value $f^{*}=0$, we have

$f(x^{k+1})-f^{*}=(1-(1-\gamma^{2}))(f(x^{k})-f^{*})$ (31)

Our algorithm (inexact optimization): We now consider the inexact active set
algorithm (denoted by ASGD). Suppose, we start with initial active set
$\mathcal{W}^{0}=\\{1\\}$ and $x^{0}=(1,1,\ldots,1)$. We take a gradient step
over $\mathcal{W}^{0}$ to obtain $x^{1}=(0,1,\ldots,1)$. Starting with
$\mathcal{W}^{k-1}$ and $x^{k}$, we randomly select $i_{k}$ from
$[p]\backslash\mathcal{W}^{k-1}$, and augment $\mathcal{W}^{k-1}$ with
$i_{k}$, i.e. $\mathcal{W}^{k}=\mathcal{W}^{k-1}\cup\\{i_{k}\\}$. Then, we
take a gradient step from $x^{k}$ over $\mathcal{W}^{k}$ to obtain $x^{k+1}$,
it is easy to see that

$x^{k+1}_{i}=\gamma x^{k}_{i},\forall i\in\mathcal{W}^{k},\quad\text{and}\quad
x^{k+1}_{i}=x^{k}_{i},\forall i\notin\mathcal{W}^{k}.$

Moreover, by induction, we have that $x^{k+1}_{1}=0$, $x^{k+1}$ contains
$(p-k-1)$ coordinates of 1, and $k$ coordinates that are
$\gamma,\gamma^{2},\ldots,\gamma^{k}$. Therefore,

$f(x^{k})=(p-k-1)+\sum_{j=1}^{k}\gamma^{2k}=(p-k-1)+\frac{\gamma^{2}-\gamma^{2k+2}}{1-\gamma^{2}}.$

Since $f^{*}=0$, we have for $k\geq 1$

$\frac{f(x^{k})-f(x^{k+1})}{f(x^{k})}=\frac{1-\gamma^{2k+2}}{p-k-1+\frac{\gamma^{2}-\gamma^{2k+2}}{1-\gamma^{2}}}\geq\frac{1-\gamma^{2}}{p-1}.$

The above inequality is tight, and it indicates

$f(x^{k+1})-f^{*}\leq\left(1-\frac{1-\gamma^{2}}{p-1}\right)(f(x^{k})-f^{*})$
(32)

Compared to the rate of gradient descent in (31), the above rate (32) has an
additional $O(\frac{1}{p})$ factor. This is similar to having the factor
$\alpha_{\\{\ell\\}}$ for the augmentation rules in Theorem 1.

##  Appendix C Additional Experiment Details

###  C.1 Real dataset details

RD1: This dataset is taken from [26] and is available from
http://ampd.epa.gov/ampd/; and has been recently used in [25, 6]. We apply the
log-transformation on all features and the response, and consider a training
set with $n=30,000$ and $d=4$. We perturb the responses with iid Gaussian
noise with mean zero and variance $\text{Var}(y)/3$.
RD2: This is taken from [32] and is the ex1029 dataset available from R
package Sleuth2. RD2 has been used in earlier work [16, 3]. We first winsorize
the data by excluding the samples that have response with score $\geq 2$; and
consider a training set with $n=15,000$ and $d=4$. Following [3], we apply the
transformation $x\mapsto 1.2^{x}$ to the education variable.

###  C.2 Algorithm Parameters

The tolerance for violations of constraints is set as $10^{-4}$ ($10^{-8}$ for
the second stage if any). For inexact optimization, the tolerance for the
relative objective change is $10^{-6}$ and the maximum number of PGD
iterations is taken to be 5. For exact optimization, the violation tolerance
is $10^{-7}$ and the maximum number of PGD iterations for the sub-problems is
$3,000$.

In Rule 1, we take $P=1$; for Rule 2, we take $K=n$; in Rule 4, we take $M=4n$
and $K=n$.