Curvature-Exploiting Acceleration of Elastic Net Computations

TL;DR

This paper presents a second-order optimization method for the elastic net that leverages curvature information to achieve faster convergence, especially on ill-conditioned datasets, with strong theoretical guarantees and practical benefits.

Contribution

It introduces a novel curvature-exploiting technique for elastic net optimization, combining second-order information with variance reduction and momentum acceleration.

Findings

Achieves faster run times than first-order methods.

Demonstrates significant practical improvements on ill-conditioned data.

Provides theoretical guarantees for the proposed method.

Abstract

This paper introduces an efficient second-order method for solving the elastic net problem. Its key innovation is a computationally efficient technique for injecting curvature information in the optimization process which admits a strong theoretical performance guarantee. In particular, we show improved run time over popular first-order methods and quantify the speed-up in terms of statistical measures of the data matrix. The improved time complexity is the result of an extensive exploitation of the problem structure and a careful combination of second-order information, variance reduction techniques, and momentum acceleration. Beside theoretical speed-up, experimental results demonstrate great practical performance benefits of curvature information, especially for ill-conditioned data sets.

Figures

Tables

Table 1. Table 1: Summary of different algorithms solving the elastic net problem. Here, κ 𝜅 \kappa and κ ~ ~ 𝜅 \tilde{\kappa} are conventional condition numbers satisfying κ ≤ κ ~ 𝜅 ~ 𝜅 \kappa\leq\tilde{\kappa} , while κ 𝑯 subscript 𝜅 𝑯 \kappa_{{}_{\bm{H}}} is a new condition number defined w.r.t the 𝑯 𝑯 \bm{H} -norm. When 𝑯 𝑯 \bm{H} is an approximate Hessian of the ridge loss, κ 𝑯 subscript 𝜅 𝑯 \kappa_{{}_{\bm{H}}} is often much smaller than κ ~ ~ 𝜅 \tilde{\kappa} , especially on practical data sets.

Algorithm	Time complexity	2nd-order \bigstrut
PGD	$O\left(dn\kappa \log \frac{1}{ϵ}\right)$	no
FISTA	$O\left(dn\sqrt{\kappa }\log \frac{1}{ϵ}\right)$	no
ProxSVRG	$O\left(d\left(n+\tilde{\kappa }\right)\log \frac{1}{ϵ}\right)$	no
Katyusha	$O\left(d\left(n+\sqrt{n\tilde{\kappa }}\right)\log \frac{1}{ϵ}\right)$	no
Ours	$O\left(d\left(n+{\kappa }_{{}_{𝑯}}\right)\log \frac{1}{ϵ}\right)$	yes

Table 2. Table 2: Brief summary of data sets and parameters

Data set	$d$	$n$	${\gamma }_1$	$r$ \bigstrut
gisette-scale	5,000	6,000	${10}^{-3}$	40
australian	14	690	${10}^{-3}$	5
cina0	132	16,033	${10}^{-4}$	20
realsim	20,958	72,309	${10}^{-4}$	50

Equations

$$\displaystyle\underset{\bm{x}\in\mathbbmss{R}^{d}}{\mbox{minimize}}\,\bigg{\{}\frac{1}{2n}\left\|{\bm{A}\bm{x}-\bm{b}}\right\|_{2}^{2}+\frac{\gamma_{2}}{2}\left\|{\bm{x}}\right\|_{2}^{2}+\gamma_{1}\left\|{\bm{x}}\right\|_{1}\bigg{\}},$$

$$\underset{\mathbf{x}\in \text{mathbbmss}{R}^d}{\min }F\left(\mathbf{x}\right)\triangleq f\left(\mathbf{x}\right)+h\left(\mathbf{x}\right),$$

$${\Vert \nabla f_i\left(\mathbf{x}\right)-\nabla f_i\left(\mathbf{y}\right)\Vert }_{{\mathbf{H}}^{-1}}\le L_i{\Vert \mathbf{x}-\mathbf{y}\Vert }_{\mathbf{H}},\forall \mathbf{x},\mathbf{y}\in \text{mathbbmss}{R}^d.$$

$${\Vert \nabla f\left(\mathbf{x}\right)-\nabla f\left(\mathbf{y}\right)\Vert }_{{\mathbf{H}}^{-1}}\le L{\Vert \mathbf{x}-\mathbf{y}\Vert }_{\mathbf{H}}$$

$$f\left(\mathbf{y}\right)\le f\left(\mathbf{x}\right)+⟨\nabla f\left(\mathbf{x}\right),\mathbf{y}-\mathbf{x}⟩+\frac{L_{\mathrm{avg}}}{2}{\Vert \mathbf{y}-\mathbf{x}\Vert }_{\mathbf{H}}^2.$$

$$f\left(\lambda \mathbf{x}+\left(1-\lambda \right)\mathbf{y}\right)$$

$$f\left(\mathbf{y}\right)\ge f\left(\mathbf{x}\right)+⟨\nabla f\left(\mathbf{x}\right),\mathbf{y}-\mathbf{x}⟩+\frac{\mu }{2}{\Vert \mathbf{y}-\mathbf{x}\Vert }_{\mathbf{H}}^2,$$

$$\displaystyle{\operatorname{prox}}_{h}^{\bm{H}}\left(\bm{y}\right)=\operatorname*{argmin}_{\bm{x}\in\mathbbmss{R}^{d}}\big{\{}h\left(\bm{x}\right)+\frac{1}{2}\left\|{\bm{x}-\bm{y}}\right\|_{\bm{H}}^{2}\big{\}}.$$

$$\displaystyle{\operatorname{prox}}_{h}\left(\bm{y}\right)=\operatorname*{argmin}_{\bm{x}\in\mathbbmss{R}^{d}}\big{\{}h\left(\bm{x}\right)+\frac{1}{2}\left\|{\bm{x}-\bm{y}}\right\|_{2}^{2}\big{\}}.$$

$$\mathbf{H}\left(\mathbf{x}-{\mathrm{prox}}_h^{\mathbf{H}}\left(\mathbf{x}\right)\right)\in \partial h\left({\mathrm{prox}}_h^{\mathbf{H}}\left(\mathbf{x}\right)\right).$$

$${\Vert {\mathrm{prox}}_h^{\mathbf{H}}(\mathbf{x})-{\mathrm{prox}}_h^{\mathbf{H}}(\mathbf{y})\Vert }_{\mathbf{H}}\le {\Vert \mathbf{x}-\mathbf{y}\Vert }_{\mathbf{H}}\forall \mathbf{x},\mathbf{y}\in \mathrm{dom}h.$$

$$\displaystyle h\left(\bm{x}^{+}\right)+\frac{1}{2\eta}\left\|{\bm{x}^{+}-\bm{y}}\right\|_{\bm{H}}^{2}\leq\underset{\bm{x}\in\mathbbmss{R}^{d}}{\min}\big{\{}h\left(\bm{x}\right)+\frac{1}{2\eta}\left\|{\bm{x}-\bm{y}}\right\|_{\bm{H}}^{2}\big{\}}+\epsilon.$$

$$\Vert \mathbf{A}-{\mathbf{A}}_r{\Vert }_2$$

$$\displaystyle\big{|}{\bm{\bm{u}}}_{i}^{\top}\bm{A}\bm{A}^{\top}{\bm{\bm{u}}}_{i}-\bar{\bm{u}}_{i}^{\top}\bm{A}\bm{A}^{\top}\bar{\bm{u}}_{i}\big{|}$$

$$\mathbf{H}={\mathbf{V}}_r\left({\mathbf{\Sigma }}_r^2+{\gamma }_2\mathbf{I}\right){\mathbf{V}}_r^{\top }+\left({\sigma }_r^2+{\gamma }_2\right)\left(\mathbf{I}-{\mathbf{V}}_r{\mathbf{V}}_r^{\top }\right).$$

$${\mathbf{H}}^{-1}={\mathbf{V}}_r{\left({\mathbf{\Sigma }}_r^2+{\gamma }_2\mathbf{I}\right)}^{-1}{\mathbf{V}}_r^{\top }+\frac{1}{{\sigma }_r^2+{\gamma }_2}\left(\mathbf{I}-{\mathbf{V}}_r{\mathbf{V}}_r^{\top }\right),$$

$${\kappa }_{{}_{\mathbf{H}}}=\frac{L_{\mathrm{avg}}}{\mu }=\frac{\frac{1}{n}{\sum }_{i=1}^n{L}_i}{\mu }.$$

$$\frac{1}{n}\sum _{i=1}^n\underset{f_i(\mathbf{x})}{\underbrace{\frac{1}{2}{\left({\mathbf{a}}_i^{\top }\mathbf{x}-{\mathbf{b}}_i\right)}^2+{\gamma }_2{\Vert \mathbf{x}\Vert }_2^2}}.$$

$$L_i$$

$$\mu$$

$$L_i=\Vert {\mathbf{a}}_i{\mathbf{a}}_i^{\top }+{\gamma }_2\mathbf{I}{\Vert }_2^2\text{mbox}and\mu ={\lambda }_d\left(\mathbf{C}+{\gamma }_2\mathbf{I}\right).$$

$${\mathsf{d}}_{\lambda }=\sum _{i=1}^d\frac{{\lambda }_i}{{\lambda }_i+\lambda }.$$

$${\lambda }_1\left({\mathbf{H}}^{-1/2}\left(\mathbf{C}+{\gamma }_2\mathbf{I}\right){\mathbf{H}}^{-1/2}\right)\le 17$$

$$\frac{{\gamma }_2}{19\left({\lambda }_r+{\gamma }_2\right)}\le {\lambda }_d\left({\mathbf{H}}^{-1/2}\left(\mathbf{C}+{\gamma }_2\mathbf{I}\right){\mathbf{H}}^{-1/2}\right)\le 2.$$

$${\kappa }_{{}_{\mathbf{H}}}\le \min \left(\frac{{\mathsf{d}}_{{\gamma }_2}}{{\gamma }_2},\frac{r{\lambda }_r+{\sum }_{i>r}{\lambda }_i}{{\gamma }_2}+d\right).$$

$$\frac{{\sum }_{i\le r}{\lambda }_i+{\sum }_{i>r}{\lambda }_i}{r{\lambda }_r+{\sum }_{i>r}{\lambda }_i},$$

$${\mathbf{v}}_k=\left(\nabla f_{i_k}\left({\mathbf{x}}_k\right)-\nabla f_{i_k}\left({\tilde{\mathbf{x}}}_s\right)\right)/\left(n{p}_{i_k}\right)+\nabla f\left({\tilde{\mathbf{x}}}_s\right)$$

$${\mathbf{x}}_{k+1}={\mathrm{prox}}_{\eta h}\left({\mathbf{x}}_k-\eta {\mathbf{v}}_k\right),$$

$${\mathbf{x}}_{k+1}={\mathrm{prox}}_{\eta h}^{\mathbf{H}}\left({\mathbf{x}}_k-\eta {\mathbf{H}}^{-1}{\mathbf{v}}_k\right),$$

$$\displaystyle\bm{v}_{k}=\frac{1}{b}\sum_{i_{k}\in\mathcal{B}_{k}}\big{(}\nabla f_{i_{k}}\left(\bm{x}_{k}\right)-\nabla f_{i_{k}}\left(\tilde{\bm{x}}_{s}\right)\big{)}/\left(np_{i_{k}}\right)+\nabla f\left(\tilde{\bm{x}}_{s}\right),$$

$$\underset{\mathbf{x}\in \text{mathbbmss}{R}^d}{\text{mbox}minimize}\frac{1}{2\eta }{\Vert \mathbf{x}-{\mathbf{y}}_k+\eta {\mathbf{H}}^{-1}{\mathbf{v}}_k\Vert }_{\mathbf{H}}^2+h\left(\mathbf{x}\right).$$

$$V_k=F\left({\mathbf{x}}_k\right)-F\left({\mathbf{x}}^{\star }\right)+\frac{\mu }{2}{\Vert {\mathbf{z}}_k-{\mathbf{x}}^{\star }\Vert }_{\mathbf{H}}^2.$$

Full text

1 Introduction

Lasso, ridge and elastic net regression are fundamental problems in statistics and machine learning, with countless applications in science and engineering [33]. Elastic net regression amounts to solving the following convex optimization problem

[TABLE]

for given data matrices $\bm{A}\in\mathbbmss{R}^{n\times d}$ and $\bm{b}\in\mathbbmss{R}^{n}$ and regularization parameters $\gamma_{1}$ and $\gamma_{2}$. Setting $\gamma_{1}=0$ results in ridge regression, $\gamma_{2}=0$ yields lasso and letting $\gamma_{1}=\gamma_{2}=0$ reduces the problem to the classical least-squares. Lasso promotes sparsity of the optimal solution, which sometimes helps to improve interpretability of the results. Adding the additional $l_{2}$-regularizer helps to improve the performance when features are highly correlated [28, 33].

The convergence rates of iterative methods for solving (1) are typically governed by the condition number of the Hessian matrix of the ridge loss, $\bm{C}+\gamma_{2}\mathbf{I}$, where $\bm{C}=\frac{1}{n}\bm{A}^{\top}\bm{A}$ is the sample correlation matrix. Real-world data sets often have few dominant features, while the other features are highly correlated with the stronger ones [28, 13]. This translates to a rapidly decaying spectrum of $\bm{C}$. In this paper, we demonstrate how this property can be exploited to reduce the effect of ill-conditioning and to design faster algorithms for solving the elastic net regression problem (1).

1.1 Related work

Over the past few years, there has been a great attention to developing efficient optimization algorithms for minimizing composite objective functions

[TABLE]

where $f\left(\bm{x}\right)=\frac{1}{n}\sum_{i=1}^{n}f_{i}\left(\bm{x}\right)$ is a finite sum of smooth and convex component functions $f_{i}\left(\bm{x}\right)$, and $h\left(\bm{x}\right)$ is a possibly non-smooth convex regularizer. In machine learning applications, the function $f$ typically models the empirical data loss and the regularizer $h$ is used to promote desired properties of a solution. For example, the elastic net objective can fit to this form with $f_{i}(\bm{x})=\frac{1}{2}\left(\bm{a}_{i}^{\top}\bm{x}-\bm{b}_{i}\right)^{2}+\gamma_{2}\left\|{\bm{x}}\right\|_{2}^{2},$ and $h(\bm{x})=\gamma_{1}\left\|{\bm{x}}\right\|_{1}.$

1.1.1 First-Order methods

Standard deterministic first-order methods for solving (2), such as proximal gradient descent, enjoy linear convergence for strongly convex objective functions and are able to find an $\epsilon$-approximate solution in time $O\left(dn\kappa\log\frac{1}{\epsilon}\right)$, where $\kappa$ is the condition number of $f$. This runtime can be improved to $O\left(dn\sqrt{\kappa}\log\frac{1}{\epsilon}\right)$ if it is combined with Nesterov acceleration [4, 23]. However, the main drawback of these methods is that they need to access the whole data set in every iteration, which is too costly in many machine learning tasks.

For large-scale problems, methods based on stochastic gradients have become the standard choice for solving (2). Many linearly convergent proximal methods such as, SAGA [9] and Prox-SVRG [31], have been introduced and shown to outperform standard first-order methods under certain regularity assumptions. These methods improve the time complexity to $O\left(d\left(n+\tilde{\kappa}\right)\log\frac{1}{\epsilon}\right)$, where $\tilde{\kappa}$ is a condition number satisfying $\tilde{\kappa}\geq\kappa$. When the component functions do not vary substantially in smoothness, $\tilde{\kappa}\approx\kappa$, and this complexity is far better than those of deterministic methods above. By exploiting Nesterov momentum in different ways (see, e.g., [11, 19, 2, 8]), one can improve the complexity to O\big{(}d\big{(}n+\sqrt{n\tilde{\kappa}}\big{)}\log\frac{1}{\epsilon}\big{)}, which is also optimal for this class of problems [30].

1.1.2 Second-order methods

Second-order methods are known to have superior performance compared to their first-order counterparts both in theory and practice, especially when the problem at hand is highly nonlinear and/or ill-conditioned. However, such methods often have very high computational cost per iteration. Recently, there has been an intense effort to develop algorithms which use second-order information with a more reasonable computational burden (see, e.g., [10, 21, 26, 1, 31, 32, 6] and references therein). Those methods use techniques such as random sketching, matrix sampling, and iterative estimation to construct an approximate Hessian matrix. Local and global convergence guarantees have been derived under various assumptions. Although many experimental results have shown excellent performance of those methods on many machine learning tasks, current second-order methods for finite-sum optimization tend to have much higher time-complexities than their first-order counterparts (see [32] for a detailed comparison).

Apart from having high time complexities, none of the methods cited above have any guarantees in the composite setting since their analyses hinge on differentiability of the objective function. Instead, one has to rely on methods that build on proximal Newton updates (see, e.g., [18, 20, 12, 25]). However, these methods still inherit the high update and storage costs of conventional second-order methods or require elaborate tuning of several parameters and stopping criteria depending on a phase transition which occurs in the algorithm.

1.1.3 Ridge regression

For the smooth ridge regression problem, the authors in [13] have developed a preconditioning method based on linear sketching which, when coupled with SVRG, yields a significant speed-up over stochastic first-order methods. This is a rare second-order method that has a comparable or even better time complexity than stochastic first-order methods. More precisely, it has a guaranteed running time of $O(d(n+\kappa_{{}_{\bm{H}}})\log\frac{1}{\epsilon})$, where $\kappa_{{}_{\bm{H}}}$ is a new condition number that can be dramatically smaller than $\tilde{\kappa}$, especially when the spectrum of $\bm{C}$ decays rapidly. When $d\ll n$, the authors in [29] combine sub-sampled Newton methods with the mini-batch SVRG to obtain some further improvements.

1.2 Contributions

Recently, the work [3] shows that under some mild algorithmic assumptions, and if the dimension is sufficiently large, the iteration complexity of second-order methods for smooth finite-sum problems composed of quadratics is no better than first-order methods. Therefore, it is natural to ask whether one can develop a second-order method for solving the elastic net problem which has improved practical performance but still enjoys a strong worst-case time complexity like the stochastic first-order methods do? It should be emphasized that due to the non-smooth objective, achieving this goal is much more challenging than for ridge regression. The preconditioning approach in [13] is not applicable, and the current theoretical results for second-order methods are not likely to offer the desired running time.

In this paper, we provide a positive answer to this question. Our main contribution is the design and analysis of a simple second-order method for the elastic net problem which has a strong theoretical time complexity and superior practical performance. The convergence bound adapts to the problem structure and is governed by the spectrum and a statistical measure of the data matrix. These quantities often yield significantly stronger time complexity guarantees for practical datasets than those of stochastic first-order methods (see Table 1). To achieve this, we first leverage recent advances in randomized low-rank approximation to generate a simple, one-shot approximation of the Hessian matrix. We then exploit the composite and finite-sum structure of the problem to develop a variant of the ProxSVRG method that builds upon Nesterov’s momentum acceleration and inexact computations of scaled proximal operators, which may be of independent interest. We provide a simple convergence proof based on an explicit Lyapunov function, thus avoiding the use of sophisticated stochastic estimate sequences.

2 Preliminaries and Notation

Vectors are indicated by bold lower-case letters, and matrices are denoted by bold upper-case letters. We denote the dot product between $\bm{x}$ and $\bm{y}$ by $\left\langle\bm{x},\bm{y}\right\rangle=\bm{x}^{\top}\bm{y}$, and the Euclidean norm of $\bm{x}$ by $\left\|{\bm{x}}\right\|_{2}=\sqrt{\left\langle\bm{x},\bm{x}\right\rangle}$. For a symmetric positive definite matrix $\bm{H}$, $\left\langle\bm{x},\bm{y}\right\rangle_{\bm{H}}=\bm{x}^{\top}\bm{H}\bm{y}$ is the $\bm{H}$-inner product of two vectors $\bm{x}$ and $\bm{y}$ and $\left\|{\bm{x}}\right\|_{\bm{H}}=\sqrt{\left\langle\bm{x},\bm{x}\right\rangle_{\bm{H}}}$ is the Mahalanobis norm of $\bm{x}$. We denote by $\lambda_{i}\left(\bm{A}\right)$ the $i$th largest eigenvalue of $\bm{A}$. Finally, $\lambda_{i}$ denotes the $i$th largest eigenvalue of the correlation matrix $\bm{C}$.

In the paper, we shall frequently use the notions of strong convexity and smoothness in the $\bm{H}$-norm, introduced in the next two assumptions.

Assumption 1**.**

The function $h\left(\bm{x}\right)$ is lower-semicontinous and convex and $\mathrm{dom}\,h:=\{\bm{x}\in\mathbbmss{R}^{d}\,|\,h\left(\bm{x}\right)<\infty\}$, is closed. Each function $f_{i}$ is $L_{i}$-smooth w.r.t the $\bm{H}$-norm, i.e, there exists a positive constant $L_{i}$ such that

[TABLE]

Assumption 1 implies that $\nabla f$ is $L$-Lipschitz:

[TABLE]

for some $L\leq L_{\mathrm{avg}}=\frac{1}{n}\sum_{i=1}^{n}L_{i}$. As a consequence, we have the following bound:

[TABLE]

Assumption 2**.**

The function $f\left(\bm{x}\right)$ is $\mu$-strongly convex w.r.t the $\bm{H}$-norm, i.e, there exists a positive constant $\mu$ such that

[TABLE]

holds for all $x,y\in\mathbbmss{R}^{n}$ and $\lambda\in\left[0,1\right]$.

Assumption 2 is equivalent to the requirement that

[TABLE]

holds for all $\bm{x},\bm{y}\in\mathbbmss{R}^{n}$. We will use both of these definitions of strong convexity in our proofs.

At the core of our method is the concept of scaled proximal mappings, defined as follows:

Definition 1** (Scaled Proximal Mapping).**

For a convex function $h$ and a symmetric positive definite matrix $\bm{H}$, the scaled proximal mapping of $h$ at $\bm{x}$ is

[TABLE]

The scaled proximal mappings generalize the conventional ones:

[TABLE]

However, while many conventional prox-mappings admit analytical solutions, this is almost never the case for scaled proximal mappings. This makes it hard to extend efficient first-order proximal methods to second-order ones. Fortunately, scaled proximal mappings do share some key properties with the conventional ones. We collect a few of them in the following result:

Property 1** ([18]).**

The following properties hold:

1. ${\operatorname{prox}}_{h}^{\bm{H}}\left(\bm{x}\right)$ exists and is unique for $\bm{x}\in\mathrm{dom}\,h$.

2. Let $\partial h\left(\bm{x}\right)$ be the subdifferential of $h$ at $\bm{x}$, then

[TABLE]

3. ${\operatorname{prox}}_{h}^{\bm{H}}\left(\cdot\right)$ is non-expansive in the $\bm{H}$-norm:

[TABLE]

Finally, in our algorithm, it will be enough to solve (3) approximately in the following sense:

Definition 2** (Inexact subproblem solutions).**

We say that $\bm{x}^{+}\in\mathbbmss{R}^{d}$ is an $\epsilon$-optimal solution to (3) if

[TABLE]

3 Building Block 1: Randomized Low-Rank Approximation

The computational cost of many Newton-type methods is dominated by the time required to compute the update direction $\bm{d}=\bm{H}^{-1}\bm{g}$ for some vector $\bm{g}\in\mathbbmss{R}^{d}$ and approximate Hessian $\bm{H}$. A naive implementation using SVD would take $O\left(nd^{2}\right)$ flops, which is prohibitive for large-scale data sets. A natural way to reduce this cost is to use truncated SVD. However, standard deterministic methods such as the power method and the Lanzcos method have run times that scale inversely with the gap between the eigenvalues of the input matrix. This gap can be arbitrarily small for practical data sets, thereby preventing us from obtaining the desired time complexity. In contrast, randomized sketching schemes usually admit gap-free run times [15]. However, unlike other methods, the block Lanczos method, detailed in Algorithm 1, admits both fast run times and strong guarantees on the errors between the true and the computed approximate singular vectors. This property turns out to be critical for deriving bounds on the condition number of the elastic net.

Proposition 1** ([22]).**

Assume that ${\bm{U}}_{r}$, ${\bm{\Sigma}}_{r}$, and ${\bm{V}}_{r}$ are matrices generated by Algorithm 1. Let ${\bm{A}}_{r}={\bm{U}}_{r}{\bm{\Sigma}}_{r}{\bm{V}}_{r}^{\top}=\sum_{i=1}^{r}{\sigma}_{i}{\bm{\bm{u}}}_{i}{\bm{\bm{v}}}_{i}^{\top}$ and let $\bm{A}=\sum_{i=1}^{d}\bar{\sigma}_{i}\bar{\bm{u}}_{i}\bar{\bm{v}}_{i}^{\top}$ be the SVD of $\bm{A}$. Then, the following bounds hold with probability at least $9/10$:

[TABLE]

The total running time is $O\left(ndr\log d(\epsilon^{\prime})^{-1/2}\right)$.

Note that we only run Algorithm 1 once and $\epsilon^{\prime}=1/2$ is sufficient in our work. Thus, the computational cost of this step is negligible, in theory and in practice.

3.1 Aproximating the Hessian

In this work, we consider the following approximate Hessian matrix of the ridge loss:

[TABLE]

Here, the first term is a natural rank $r$ approximation of the true Hessian, while the second term is used to capture information in the subspace orthogonal to the column space of $\bm{V}_{r}$. The inverse of $\bm{H}$ in (6) admits the explicit expression

[TABLE]

so the evaluation of $\bm{H}^{-1}\bm{x}$ has time complexity $O\left(rd\right)$.

3.2 Bounding the Condition Number

We now turn our attention to studying how the approximate Hessian affects the relevant condition number of the elastic net problem. We first introduce a condition number that usually determines the iteration complexity of stochastic first-order methods under non-uniform sampling.

Definition 3**.**

The average condition number of (1.1) is

[TABLE]

For the elastic net problem (1), the smooth part of the objective is the ridge loss

[TABLE]

Since we define smoothness and strong convexity of $f_{i}(\bm{x})$ in the $\bm{H}$-norm, the relevant constants are

[TABLE]

For comparison, we also define the conventional condition number $\tilde{\kappa}$, which characterizes the smoothness and strong convexity of $f_{i}(\bm{x})$ in the Euclidean norm. In this case $\tilde{\kappa}=\sum_{i}L_{i}/(n\mu)$, where

[TABLE]

It will become apparent that $\kappa_{{}_{\bm{H}}}$ can be expressed in terms of a statistical measure of the ridge loss and that it may be significantly smaller than $\tilde{\kappa}$. We start by introducing a statistical measure that has been widely used in the analysis of ridge regression (see, e.g., [16] and the references therein).

Definition 4** (Effective Dimension).**

For a positive constant $\lambda$, the effective dimension of $\bm{C}$ is defined as

[TABLE]

The effective dimension generalizes the ordinary dimension and satisfies $\mathsf{d}_{\lambda}\leq d$ with equality if and only if $\lambda=0$. It is typical that when $\bm{C}$ has a rapidly decaying spectrum, most of the $\lambda_{i}$’s are dominated by $\lambda$, and hence $\mathsf{d}_{\lambda}$ can be much smaller than $d$.

The following lemma bounds the eigenvalues of the matrix $\bm{H}^{-1/2}\left(\bm{C}+\gamma_{2}\mathbf{I}\right)\bm{H}^{-1/2}$, which can be seen as the effective Hessian matrix.

Lemma 1** (​​[13]).**

Invoking Algorithm 1 with data matrix $\frac{1}{\sqrt{n}}\bm{A}$, target rank $r$, and target precision $\epsilon^{\prime}=1/2$, it holds with probability at least $9/10$ that

[TABLE]

Equipped with Lemma 1, we can now connect $\kappa_{{}_{\bm{H}}}$ with $\mathsf{d}_{\lambda}$ using the following result.

Theorem 1**.**

With probability at least $9/10$, the following bound holds up to a multiplicative constant:

[TABLE]

Proof.

See Appendix A. ∎

Since $\tilde{\kappa}=(\sum_{i}\lambda_{i}+d\gamma_{2})/\gamma_{2}$, $\kappa_{{}_{\bm{H}}}$ is reduced by a factor

[TABLE]

compared to $\tilde{\kappa}$. If the spectrum of $\bm{C}$ decays rapidly, then the terms $\sum_{i>r}\lambda_{i}$ are negligible and the ratio is approximately $\sum_{i\leq r}\lambda_{i}/(r\lambda_{r})$. If the first eigenvalues are much larger than $\lambda_{r}$, this ratio will be large. For example, for the australian data set [7], this ratio can be as large as $1.34\times 10^{4}$ and $1.6\times 10^{5}$ for $r=3$ and $r=4$, respectively. This indicates that it is possible to improve the iteration complexity of stochastic first-order methods if one can capitalize on the notions of strong convexity and smoothness w.r.t the $\bm{H}$-norm in the optimization algorithm. Of course, this is only meaningful if there is an efficient way to inject curvature information into the optimization process without significantly increasing the computational cost. In the smooth case, i.e., $\gamma_{1}=0$, this task can be done by a preconditioning step [13]. However, this approach is not applicable for the elastic net, and we need to make use of another building block.

4 Building Block 2: Inexact Accelerated Scaled Proximal SVRG

In this section, we introduce an inexact scaled accelerated ProxSVRG algorithm for solving the generic finite-sum minimization problem in (2). We then characterize the convergence rate of the proposed algorithm.

4.1 Description of the Algorithm

To motivate our algorithm, we first recall the ProxSVRG method from [31]: For the $s$th outer iteration with the corresponding outer iterate $\tilde{\bm{x}}_{s}$, let $\bm{x}_{0}=\tilde{\bm{x}}_{s}$ and for $k=0,2,\ldots,T-1$ do

[TABLE]

where $i_{k}$ is drawn randomly from $\{1,\ldots,n\}$ with probability $p_{i_{k}}=L_{i_{k}}/(nL_{\mathrm{avg}})$. Since we are provided with an approximate Hessian matrix $\bm{H}$, it is natural to use the following update:

[TABLE]

which can be seen as a proximal Newton step with the full gradient vector replaced by $\bm{v}_{k}$. Note that when $h(\cdot)$ is the $\ell_{1}$-penalty, ProxSVRG can evaluate (8) in time $O(d)$, while evaluating (9) amounts to solving an optimization problem. It is thus is critical to keep the number of such evaluations small, which then translates into making a sufficient progress at each iteration. A natural way to achieve this goal is to reduce the variance of the noisy gradient $\bm{v}_{k}$. This suggests to use large mini-batches, i.e., instead of using a single component function $f_{i_{k}}$, we use multiple ones to form:

[TABLE]

where $\mathcal{B}_{k}\subset\{1,\ldots,n\}$ is a set of indices with cardinality $\left|{\mathcal{B}_{k}}\right|=b$. It is easy to verify that $\bm{v}_{k}$ is an unbiased estimate of $\nabla f\left(\bm{x}_{k}\right)$. Notice that naively increasing the batch size makes the algorithm increasingly similar to its deterministic counterpart, hence inheriting a high-time complexity. This makes it hard to retain the runtime of ProxSVRG in the presence of 2nd-order information.

In the absence of second-order information and under the assumption that the proximal step is computed exactly, the work [24] introduced a method called AccProxSVRG that enjoys the same time complexity as ProxSVRG but allows for much larger mini-batch sizes. In fact, it can tolerate a mini-batch of size O\big{(}\sqrt{\tilde{\kappa}}\big{)} thanks to the use of Nesterov momentum. This indicates that an appropriate use of Nesterov momentum in our algorithm could allow for the larger mini-batches required to balance the computational cost of using scaled proximal mappings. The improved iteration complexity of the scaled proximal mappings will then give an overall acceleration in terms of wall-clock time. As discussed in [2], the momentum mechanism in AccProxSVRG fails to accelerate ProxSVRG unless $\tilde{\kappa}\geq n^{2}$. In contrast, as we will see, our algorithm will be able to accelerate the convergence also in these scenarios. In summary, our algorithm is designed to run in an inner-outer fashion as ProxSVRG with large mini-batch sizes and Nesterov momentum to compensate for the increased computational cost of subproblems. The overall procedure is summarized in Algorithm 2.

4.2 Convergence Argument

In this subsection, we will show that as long as the errors in evaluating the scaled proximal mappings are controlled in an appropriate way, the iterates generated by the outer loop of Algorithm 2 converge linearly in expectation to the optimal solution. Recall that in Step 7 of Algorithm 2, we want to find an $\epsilon_{k}$-optimal solution in the sense of (5) to the following problem:

[TABLE]

The next lemma quantifies the progress made by one inner iteration of the algorithm. Our proof builds on a Lyapunov argument using a Lyapunov function on the form:

[TABLE]

Lemma 2**.**

Let Assumptions 1–2 hold and let $\bm{x}^{\star}=\operatorname*{argmin}_{\bm{x}}F\left(\bm{x}\right)$, $\eta=1/L_{\mathrm{avg}}$ and $\tau=\sqrt{\mu/2L_{\mathrm{avg}}}$. If the mini-batch size is chosen such that $b\geq 60\sqrt{L_{\mathrm{avg}}/\mu}$, then for any $k\in\{0,\ldots,T-1\}$, there exists a vector $\bm{\xi}_{k}\in\mathbbmss{R}^{d}$ such that $\left\|{\bm{\xi}_{k}}\right\|_{\bm{H}^{-1}}\leq\sqrt{2\eta\epsilon_{k}}$ and

[TABLE]

Proof.

See Appendix C. ∎

Remark 1**.**

Our proof is direct and based on natural Lyapunov functions, thereby avoiding the use of stochastic estimate subsequences as in [24] which is already very complicated even when the subproblems are solved exactly and $\bm{H}=\mathbf{I}$. We stress that the result in Lemma 2 also holds for smaller mini-batch sizes, namely b\in\big{\{}1,\ldots,O(\sqrt{L_{\mathrm{avg}}/\mu})\big{\}}, provided that the step size $\eta$ is reduced accordingly. In favor of a simple proof, we only report the large mini-batch result here.

Equipped with Lemma 2, we can now characterize the progress made by one outer iteration of Algorithm 2.

Theorem 2**.**

Let Assumptions 1–2 hold. Suppose that the parameters $\eta$, $b$, and $\tau$ are chosen according to Lemma 2 and define $\rho=9\tau/10$. Then, if the errors in solving the subproblems satisfy

[TABLE]

for all $k\in\{0,\ldots,T-1\}$ and $T\geq(4\log c)/3\rho,$ where $c$ is a universal constant, then for every $s\in\mathbbmss{N}_{+}$,

[TABLE]

Proof.

See Appendix D. ∎

Remark 2**.**

The theorem indicates that if the errors in solving the subproblems are controlled appropriately, the outer iterates generated by Algorithm 2 converge linearly in expectation to the optimal solution. Since $V_{0}$ depends on $\bm{x}^{\star}$, it is difficult to provide a general closed-form expression for the target precisions $\epsilon_{k}$. However, we will show below that with a certain policy for selecting the initial point, it is sufficient to run the solver a constant number of iterations independently of $\bm{x}^{\star}$. We stress that the results in this section are valid for minimizing general convex composite functions (2) and not limited to the elastic net problem.

5 Warm-Start

The overall time complexity of Algorithm 2 depends strongly on our ability to solve problem (10) in a reasonable computational time. If one naively starts the solver at a random point, it may take many iterations to meet the target precision. Thus, it is necessary to have a well-designed warm-start procedure for initializing the solver. Intuitively, the current iterate $\bm{x}_{k}$ can be a reasonable starting point since the next iterate $\bm{x}_{k+1}$ should not be too far away from $\bm{x}_{k}$. However, in order to achieve a strong theoretical running time, we use a rather different scheme inspired by [19]. Let us first define the vector $\bm{u}_{k}=\bm{y}_{k}-\eta\bm{H}^{-1}\bm{v}_{k}$ for $k\in\{0,1,\ldots,T-1\}$ and the function

[TABLE]

Then, the $k$th subproblem seeks for $\bm{x}_{k+1}$ such that

[TABLE]

where $\bm{x}_{k+1}^{\star}$ is the exact solution. We consider the initialization policy

[TABLE]

which can be seen as one step of the proximal gradient method applied to $p\left(\bm{z},\bm{u}_{k-1}\right)$ starting at the current $\bm{x}_{k}$ with step size $\gamma$.

The following proposition characterizes the difference in objective realized by $\bm{z}_{0}$ and $\bm{x}_{k+1}^{\star}$.

Proposition 2**.**

Let $\bm{z}_{0}$ be defined by (14) with $\gamma=\eta/\lambda_{1}\left(\bm{H}\right)$. Let $\kappa_{\mathrm{sub}}=\lambda_{1}\left(\bm{H}\right)/\lambda_{r}\left(\bm{H}\right)$ be the condition number of the subproblems. Assume that the errors in solving the subproblems satisfy $\epsilon_{k}\leq\left(1-\rho\right)^{k}V_{0}$ for all $k\in\{0,1,\ldots,T-1\}$. Then,

[TABLE]

Proof.

See Appendix E. ∎

The proposition, together with (13), implies that it suffices to find $\bm{x}_{k+1}$ such that

[TABLE]

This is significant since one only needs to reduce the residual error by a constant factor independent of the target precision. Note also that the condition number $\kappa_{\mathrm{sub}}$ is much smaller than $\kappa\approx\lambda_{1}(\bm{H})/(\lambda_{d}+\gamma_{2})$, and computing the gradient of the smooth part of $p\left(\bm{z},\bm{u}\right)$ only takes time $O\left(rd\right)$ instead of $O\left(nd\right)$ as in the original problem. Those properties imply that the subproblems can be solved efficiently by iterative methods, where only a small (and known) constant number of iterations is needed. The next section develops the final details of convergence proof.

6 Global Time Complexity

We start with the time complexity of Algorithm 2. Let $\mathcal{T}\left(\alpha\right)$ be the number of gradient evaluations that a subproblem solver takes to reduce the residual error by a factor $\alpha$. Then, by Proposition 2, one can find an $\epsilon_{k}$-optimal solution to the $k$th subproblem by at most $\mathcal{T}\left(\kappa_{\mathrm{sub}}/\left(1-\rho\right)\right)$ gradient evaluations, where one gradient evaluation is equivalent to $O\left(d\right)$ flops. Consider the same setting of Theorem 2 and suppose that the subproblems are initialized by (14). Then, the time complexity of Algorithm 2 is given by:

[TABLE]

where the first summand is due to the full gradient evaluation at each outer loop; the second one comes from the fact that one needs $O(\sqrt{\kappa_{{}_{\bm{H}}}})$ inner iterations, each of which uses a mini-batch of size $O(\sqrt{\kappa_{{}_{\bm{H}}}})$; and the third one is the result of $O(\sqrt{\kappa_{{}_{\bm{H}}}})$ inner iterations, each of which solves a subproblem that needs $\mathcal{T}\left(\kappa_{\mathrm{sub}}/\left(1-\rho\right)\right)$ gradient evaluations. We can now put things together and state our main result.

Proposition 3**.**

Suppose that the approximate Hessian matrix $\bm{H}$ is given by (6) and that Algorithm 2 is invoked with $f\left(\bm{x}\right)=\frac{1}{2n}\left\|{\bm{A}\bm{x}-\bm{b}}\right\|_{2}^{2}+\frac{\gamma_{2}}{2}\left\|{\bm{x}}\right\|_{2}^{2}$ and $h\left(\bm{x}\right)=\left\|{\bm{x}}\right\|_{1}$. Assume further that the subproblems are solved by the accelerated proximal gradient descent method [23, 4]. Our method can find an $\epsilon$-optimal solution in time

[TABLE]

Proof.

The task reduces to evaluating the term $\mathcal{T}\left(\kappa_{\mathrm{sub}}/\left(1-\rho\right)\right)$ in (16). Recall that the iteration complexity of the accelerated proximal gradient descent method for minimizing the function $F(x)=f(x)+h(x)$, where $f$ is a smooth and strongly convex function and $h$ is a possibly non-smooth convex regularizer, initiliazed at $\bm{x}_{0}$ is given by $\sqrt{\kappa}\log\frac{F(\bm{x}_{0})-F(\bm{x}^{\star})}{\epsilon}$, where $\kappa$ is the condition number. By invoking the about result with $F(\bm{x})=p\left(\bm{x},\bm{u}_{k}\right)$, $\bm{x}^{\star}=\bm{x}_{k+1}^{\star}$, $\bm{x}_{0}=\bm{z}_{0}$, $\kappa=\kappa_{\mathrm{sub}}$, and $\epsilon$ is the right-hand side of (15), it follows that the number of iterations for each subproblem can be bounded by

[TABLE]

In addition, each iteration takes time $O\left(rd\right)$ to compute the gradient implying the time complexity

[TABLE]

Finally, since $\kappa_{{}_{\bm{H}}}\gg\kappa_{\mathrm{sub}}$, the proof is complete. ∎

We can easily recognize that this time complexity has the same form as the stochastic first-order methods discussed in Section 1.2.1 with the condition number $\tilde{\kappa}$ replaced by $\kappa_{{}_{\bm{H}}}$. It has been shown in Theorem 1 that $\kappa_{{}_{\bm{H}}}$ can be much smaller than $\tilde{\kappa}$, especially, when $\bm{C}$ has a rapidly decaying spectrum. We emphasize that the expression in (17) is available for free to us after having approximated the Hessian matrix. Hence no tuning is required to set $\epsilon_{k}$ for solving the subproblems.

7 Experimental Results

In this section, we perform numerical experiments to verify the efficacy of the proposed method on real world data sets [7, 14]. We compare our method with several well-known first-order methods: FISTA [4] with optimal step-size; Prox-SVRG [31] with epoch length $2n/b$ as suggested by the authors; Katyusha1 [2] with epoch length $2n/b$, Katyusha momentum $\tau_{2}=0.5/b$ as suggested by the author; and our method with epoch length $2n/b$. Since Katyusha1 can use a mini-batch of size $\sqrt{n}$ without slowing down the convergence, we set $b=\sqrt{n}$ for all methods except FISTA. Finally, to make a fair comparison, for each algorithm above, we tune only the step size, from the set $\eta\times\{10^{k},2\times 10^{k},5\times 10^{k}|k\in\{0,\pm 1,\pm 2\}\}$, where $\eta$ is the theoretical step size, and report the one having smallest objective value. Other hyper-parameters are set to their theory-predicted values. For Katyusha1, we also compute the largest eigenvalue of the Hessian matrix in order to set its parameter $\tau_{1}$. All methods are initialized at $\bm{0}$ and run for up to 100 epochs (passes through the full data set). For the subproblems in Algorithm 2, we just simply run FISTA with $\sqrt{\kappa_{\mathrm{sub}}}\log\kappa_{\mathrm{sub}}$ iterations as discussed in the previous section, without any further tunning steps. The value of $r$ is chosen as a small fraction of $d$ so that the preprocessing time of Algorithm 1 is negligible. Note that the available spectrum of $\bm{C}$ after running Algorithm 1 also provides an insightful way to choose $r$.

Figure 2 shows the suboptimality in objective versus the number of epochs for different algorithms solving the elastic net problem. We can see that our method systematically outperforms the others in all settings, and that there is a clear correspondence between the spectrum of $\bm{C}$ in Fig. 1 and the potential speed-up. Notably, for the australian dataset, all the first-order methods make almost no progress in the first 100 epochs, while our method can find a high-accuracy solution within tens of epochs, demonstrating a great benefit of second-order information. On the other hand, for a well-conditioned data set that does not exhibit high curvature such as real-sim, ProxSVRG is comparable to our method and even outperforms Katyusha1. This agrees with the theoretical time complexities summarized in Table 1. Finally, we demonstrate a hard instance for low-rank approximation methods via the cina0 data set, where $\bm{C}$ has a very large condition number and slowly decaying dominant eigenvalues. In this case, low-rank approximation methods with a small approximate rank $r$ may not be able to capture sufficient curvature information. We can see that even in this case curvature information helps to avoid the stagnation experienced by FISTA and ProxSVRG.

8 Conclusions

We have proposed and analyzed a novel second-order method for solving the elastic net problem. By carefully exploiting the problem structure, we demonstrated that it is possible to deal with the non-smooth objective and to efficiently inject curvature information into the optimization process without (significantly) increasing the computational cost per iteration. The combination of second-order information, fast iterative solvers, and a well-designed warm-start procedure results in a significant improvement of the total runtime complexity over popular first-order methods. An interesting direction for future research would be to go beyond the quadratic loss. We believe that the techniques developed in this work can be extended to more general settings, especially when the smooth part of the objective function is self-concordant.

Acknowledgments

This research was sponsored in part by the Knut and Alice Wallenberg Foundation and the Swedish Research Council.

Appendix A Proof of Theorem 1

First, by the definition of $L_{i}$, it holds that

[TABLE]

where the last inequality follows since

[TABLE]

The term $\bm{a}_{i}^{\top}\bm{H}^{-1}\bm{a}_{i}$ can be bounded as

[TABLE]

It can be verified that

[TABLE]

where we have used Lemma 1 in the last step. It follows that

[TABLE]

Note that

[TABLE]

which implies that

[TABLE]

On the other hand, it was shown in [13, Theorem 4] that

[TABLE]

Finally, recall that $\mu=\lambda_{d}\left(\bm{H}^{-1/2}\left(\bm{C}+\gamma_{2}\mathbf{I}\right)\bm{H}^{-1/2}\right)$. By Lemma 1, we readily have $\mu\geq c_{1}\gamma_{2}$ for some constant $c_{1}$, which concludes the proof.

Appendix B Some Useful Auxiliary Results

To facilitate the analysis, we collect some useful inequalities regarding the Mahalanobis norm that are used in the subsequent proofs.

• •

Cauchy’ inequality: $\left\langle\bm{x},\bm{y}\right\rangle\leq\left\|{\bm{x}}\right\|_{\bm{H}}\left\|{\bm{y}}\right\|_{\bm{H}^{-1}}$ for all $\bm{x},\bm{y}\in\mathbbmss{R}^{d}$

• •

Young’s inequality: $\left\langle\bm{x},\bm{y}\right\rangle\leq\frac{\alpha\left\|{\bm{x}}\right\|_{\bm{H}}^{2}}{2}+\frac{\left\|{\bm{y}}\right\|_{\bm{H}^{-1}}^{2}}{2\alpha}$ for all $\bm{x},\bm{y}\in\mathbbmss{R}^{d}$

• •

Strongly convex inequality: $\frac{\mu}{2}\left\|{\bm{x}-\bm{x}^{\star}}\right\|_{\bm{H}}^{2}\leq F\left(\bm{x}\right)-F\left(\bm{x}^{\star}\right)$, where $\bm{x}^{\star}=\operatorname*{argmin}_{\bm{x}}F\left(\bm{x}\right)$.

Appendix C Proof of Lemma 2

Before proving the lemma, we rewrite the sequences $\{\bm{x}_{k},\bm{y}_{k},\bm{z}_{k}\}_{k\geq 0}$ generated by Algorithm 2 in the following recurrence form:

[TABLE]

We also recall the following well-known three-point identity:

[TABLE]

which holds for any symmetric matrix $\bm{H}\succ\mathbf{0}$ and any vectors $\bm{a},\bm{b},\bm{c}\in\mathbbmss{R}^{d}$.

We are now ready to prove Lemma 2. From the definition of the Lyapunov function in (11), we have that

[TABLE]

where the equality $\mathrm{(a)}$ follows from (19), and $\mathrm{(b)}$ follows from(18c). Using the identity (19) again for the last term in (C), we obtain

[TABLE]

where the last equality follows from (18b). By noting that

[TABLE]

where we have used (18b) in the last equality. Then, combining (C) and (C) yields

[TABLE]

We now pay attention to the term $F\left(\bm{x}_{k+1}\right)$. Let $\bm{u}=\left(1-\tau\right)\bm{x}_{k}+\tau\bm{x}^{\star}$, then it follows from the strong convexity of $F$ that

[TABLE]

In addition, we have that

[TABLE]

where we have used (18a) and the fact that $\bm{y}_{k}-\bm{u}=\left(1-\tau\right)\left(\bm{y}_{k}-\bm{x}_{k}\right)+\tau\left(\bm{y}_{k}-\bm{x}^{\star}\right)$.

With these observations, we are now ready to bound $F\left(\bm{x}_{k+1}\right)$. In particular, by invoking Lemma 3 with $\bm{y}=\bm{y}_{t}$, $\bm{x}^{+}=\bm{x}_{t+1}$, $\bm{u}=\left(1-\tau\right)\bm{x}_{k}+\tau\bm{x}^{\star}$, $\bm{v}=\bm{v}_{t}$, $\bm{\xi}=\bm{\xi}_{k}$, $\epsilon=\epsilon_{k}$, together with (24)–(C), it follows that

[TABLE]

Substituting (C) into (C) and rearranging the terms to obtain

[TABLE]

We next bound the term $T_{1}$. By adding and subtracting the term $\bm{y}_{k}$ yields

[TABLE]

where we have used (18a) and (18b). By Young’s inequality, we have

[TABLE]

For $T_{2}$, we have

[TABLE]

Thus, combining (C)–(C) and using the fact that $\left\|{\bm{\xi}_{k}}\right\|_{\bm{H}^{-1}}\leq\sqrt{2\eta\epsilon_{k}}$ yield

[TABLE]

Invoking Lemma 4 with $\bm{u}_{k}=\left(1-\tau\right)\bm{x}_{k}+\tau\bm{x}^{\star}$ and then applying Lemma 5, we obtain

[TABLE]

Therefore, by taking the expectation on both sides of (C) and using (32), we obtain

[TABLE]

By choosing $\eta=\frac{1}{L_{\mathrm{avg}}}$, $\tau=\sqrt{\frac{\mu}{2L_{\mathrm{avg}}}}$, and $b\geq 60\sqrt{\frac{L_{\mathrm{avg}}}{\mu}}$, it is readily verified that the second and third terms on the right-hand side of (C) become nonpositive, which concludes the proof.

Appendix D Proof of Theorem 2

To begin with, let $e_{t}=\operatorname{\mathds{E}}\left\{{F\left(\bm{x}_{t}\right)-F\left(\bm{x}^{\star}\right)}\right\}$ and $\tilde{e}_{s}=\operatorname{\mathds{E}}\left\{{F\left(\tilde{\bm{x}}_{s}\right)-F\left(\bm{x}^{\star}\right)}\right\}$, then by applying inequality (2) recursively, we obtain

[TABLE]

It can be verified that

[TABLE]

Thus, from (34) and the facts that $\tau=\sqrt{\frac{\mu}{2L_{\mathrm{avg}}}}$, $\left\|{\bm{\xi}_{t}}\right\|_{\bm{H}^{-1}}\leq\sqrt{2\epsilon_{t}/L_{\mathrm{avg}}}$, it holds that

[TABLE]

where

[TABLE]

By the definition of $V_{k}$, we have $\frac{\mu}{2}\left\|{\bm{z}_{k}-\bm{x}^{\star}}\right\|_{\bm{H}}^{2}\leq V_{k}$, hence, it follows from (D) that

[TABLE]

Multiplying both sides of the above inequality by $\left(1-\tau\right)^{-k}$ yields

[TABLE]

Define $\alpha_{t}=2\sqrt{\frac{\epsilon_{t}}{\mu}}\left(1-\tau\right)^{-t/2-1}$ and $u_{t}=\left(1-\tau\right)^{-t/2}\operatorname{\mathds{E}}\left\|{\bm{z}_{t}-\bm{x}^{\star}}\right\|_{\bm{H}}$, then we can write the previous inequality as

[TABLE]

where we have separated the term $\alpha_{0}u_{0}$ from the sum and added the positive term $\alpha_{k}u_{k}$ in the last step. Note that $\alpha_{0}=\frac{2}{1-\tau}\sqrt{\frac{\epsilon_{0}}{\mu}}\leq\frac{2}{1-\tau}\sqrt{\frac{V_{0}}{\mu}}$ and $u_{0}=\left\|{\bm{x}_{0}-\bm{x}^{\star}}\right\|_{\bm{H}}\leq\sqrt{\frac{2}{\mu}V_{0}}$, hence, $u_{k}^{2}$ can be further bounded by

[TABLE]

where $c_{1}=1+\frac{\sqrt{2}}{1-\tau}$. It is readily verified that $\{B_{k}\}_{k}$ is an increasing sequence and that $\{S_{k}\triangleq\frac{2}{\mu}\left(c_{1}V_{0}+B_{k}\right)\}_{k}$ is an increasing sequence satisfying $S_{0}>u_{0}^{2}$. Therefore, by invoking Lemma 7 with $u_{t}$, $\alpha_{k}$, and $S_{k}$, we have for any $k\geq 1$ that

[TABLE]

where we have used $\sqrt{a+b}\leq\sqrt{a}+\sqrt{b}$ for any $a,b\geq 0$. Thus, for any $t\in\{0,1,\ldots,k\}$, it holds that

[TABLE]

Having upper bounds of $\operatorname{\mathds{E}}\left\|{\bm{z}_{t}-\bm{x}^{\star}}\right\|_{\bm{H}}$, we can now substitute them into (D) to get

[TABLE]

where the equality follows from the definition of $\alpha_{t}$, and in the last step, we have added a positive term $\alpha_{k}$ in the first sum as well as added and subtracted $\alpha_{0}$ in the second sum. Note that the function $\sqrt{1+x}-x$ is decreasing on the interval $\left[0,\infty\right)$, it follows that

[TABLE]

We thus have

[TABLE]

where in the last step, we have applied the inequality $\left(a+b\right)^{2}\leq\left(1+\beta\right)a^{2}+\left(1+\frac{1}{\beta}\right)b^{2}$ twice with $\beta=9$ and $\beta=1$, respectively. Using the definitions of $B_{k}$ and $\alpha_{t}$, the above inequality can be rewritten as

[TABLE]

With our choice of $\epsilon_{t}$, it follows that

[TABLE]

Thus, the last term on the right-hand side of (D) can be bounded by

[TABLE]

where the last inequality follows since the function $\sqrt{1-x}+x/2$ is decreasing on the interval $\left[0,1\right]$. Similarly, for the second term in (D), we have

[TABLE]

Thus, we can further bound $\operatorname{\mathds{E}}V_{k}$ as

[TABLE]

where $c_{2}=20+\frac{160\left(1-\rho\right)}{\left(\tau-\rho\right)^{2}\left(1-\tau\right)^{2}}+\frac{50}{9\left(\tau-\rho\right)}$ and we have also used the fact that $1-\tau<1-\rho$.

Now, if we let $\delta_{k}$ be the right-hand side of the above inequality, then it holds that

[TABLE]

where the third equality follows since by definition, $e_{k-1}\leq\operatorname{\mathds{E}}V_{k-1}\leq\delta_{k-1}$, and since $\tilde{e}_{s}=e_{0}\leq V_{0}$, and the last inequality follows since $\rho=0.9\tau$. By applying the above inequality recursively, we obtain

[TABLE]

where $c_{3}=c_{2}+2/9$. Therefore, when $T\geq\frac{4}{3\rho}\log\left(27c_{3}\right)$, it holds that

[TABLE]

Finally, using the definition of $V_{T}$, we obtain

[TABLE]

completing the proof.

Appendix E Proof of Proposition 2

Recall that

[TABLE]

which can be seen as one step of the proximal gradient method applied to $p\left(\bm{z},\bm{u}_{k-1}\right)$ starting at the current $\bm{x}_{k}$ with

[TABLE]

By the optimality condition of $\bm{z}_{0}$, we have

[TABLE]

Since $\partial p\left(\bm{z}_{0},\bm{u}_{k-1}\right)=\partial h\left(\bm{z}_{0}\right)+\frac{1}{\eta}\bm{H}\left(\bm{z}_{0}-\bm{u}_{k-1}\right)$, it follows that

[TABLE]

We see that $\bm{\zeta}_{k}$ is independent of the second argument of $p\left(\bm{z}_{0},\cdot\right)$, hence we also have $\bm{\zeta}_{k}\in\partial p\left(\bm{z}_{0},\bm{u}_{k}\right)$. Since $p\left(\bm{z},\cdot\right)$ is strongly convex in $\bm{z}$ with parameter $\lambda_{r}\left(\bm{H}\right)/\eta$ and $\bm{\zeta}_{k}\in\partial p\left(\bm{z}_{0},\bm{u}_{k}\right)$, it holds that

[TABLE]

where the last inequality follows from Young’s inequality. We next bound $\left\|{\bm{\zeta}_{k}}\right\|^{2}$ via

[TABLE]

where the last step follows from the definition of $\gamma$ and the fact that

[TABLE]

Thus, combining (E) and (E) yields

[TABLE]

Note that the quantity $\frac{1}{\gamma}\left(\bm{x}_{k}-\bm{z}_{0}\right)$ is nothing but the gradient mapping of $p\left(\bm{z}_{0},\bm{u}_{k-1}\right)$, hence by [23, Theorem 1],

[TABLE]

We thus have

[TABLE]

as desired.

Appendix F Proof of Auxiliary Lemmas

Recall that at each step of Algorithm 2, we wish to solve the following problem:

[TABLE]

The following lemma is a generalization of [17, Lemma 3] and [31, Lemma 3] to account for inexactness in the evaluation of the proximal operator, and a Mahalanobis norm.

Lemma 3**.**

Let Assumptions 1–2 hold. For any $\bm{y}\in\mathrm{dom}h$, and $\bm{v}\in\mathbbmss{R}^{d}$, let $\bm{x}^{+}$ be an $\epsilon$-optimal solution to problem (40) in the sense of (5), where $\eta\in\left(0,\frac{1}{L_{\mathrm{avg}}}\right]$ is a constant step-size. Denote $\Delta=\bm{v}-\nabla f\left(\bm{y}\right)$, then, there exists a vector $\bm{\xi}\in\mathbbmss{R}^{d}$ such that $\left\|{\bm{\xi}}\right\|_{\bm{H}^{-1}}\leq\sqrt{2\eta\epsilon}$, and it holds for any $\bm{u}\in\mathbbmss{R}^{d}$ that

[TABLE]

Proof.

Since $\bm{x}^{+}$ is an $\epsilon$-optimal solution to problem (40), by Lemma 6, there exists a vector $\bm{\xi}$, satisfying $\left\|{\bm{\xi}}\right\|_{\bm{H}^{-1}}\leq\sqrt{2\eta\epsilon}$, such that

[TABLE]

Since $h$ is convex, for any vector $\bm{\zeta}\in\partial_{\epsilon}h\left(\bm{x}^{+}\right)$, it holds for all $\bm{u}\in\mathbbmss{R}^{d}$ that

[TABLE]

Therefore, it follows from (41) that

[TABLE]

By the smoothness of $f$, we have for $\eta\in\left(0,\frac{1}{L_{\mathrm{avg}}}\right]$ that

[TABLE]

where the last equality follows from adding and subtracting the term $\bm{u}$ and expanding the norm squared. Since $f$ is $\mu$-strongly convex w.r.t the $\bm{H}$-norm, we have

[TABLE]

It follows that

[TABLE]

where the equality follows from adding and subtracting the terms $\bm{v}$ and $\bm{x}^{+}$ in the first and second inner products, respectively; and the last inequality follows from (42). This completes the proof of Lemma 3. ∎

Lemma 4**.**

Let Assumptions 1–2 hold. Let $\bm{x}_{k+1}$ be an $\epsilon_{k}$-optimal solution to the subproblem in Step 7 of Algorithm 2. Let $\bm{u}_{k}$ be any vector in $\mathbbmss{R}^{d}$ that is independent of the mini-batch $\mathcal{B}_{k}$, then, it holds for any $\eta\in\left(0,\frac{1}{L_{\mathrm{avg}}}\right]$ that

[TABLE]

Proof.

Let $\bm{x}_{k+1}^{\star}$ be the exact solution to the subproblem in Step 7 of Algorithm 2, i.e, $\epsilon_{k}=0$. Since the objective function $q$ defined in (40) is $\frac{1}{\eta}$-strongly convex w.r.t the $\bm{H}$-norm, we have

[TABLE]

Therefore, one can write $\bm{x}_{k+1}$ as $\bm{x}_{k+1}=\bm{x}_{k+1}^{\star}+\bm{\zeta}_{k}$ for some vector $\bm{\zeta}_{k}\in\mathbbmss{R}^{d}$ satisfying $\left\|{\bm{\zeta}_{k}}\right\|_{\bm{H}}\leq\sqrt{2\eta\epsilon_{k}}$. If we define the following virtual iterate:

[TABLE]

then

[TABLE]

where the last step follows from Cauchy’s inequality. Note that

[TABLE]

where the first inequality follows from the triangle inequality and the last one follows from Property 1. Therefore, it holds that

[TABLE]

Denote the filtration $\mathcal{F}_{k}$ by $\mathcal{F}_{k}=\{\bm{f}_{0},\bm{f}_{1},\ldots,\bm{f}_{k}\}$, where $\bm{f}_{i}$ denotes all the randomness incurring at time $i$ for all $i\in\{1,\ldots,k\}$. Note that the triple $\left(\bm{x}_{k},\bm{y}_{k},\bm{z}_{k}\right)$ depends on $\mathcal{F}_{k-1}$, but not on $\bm{f}_{k}$, we thus have

[TABLE]

Taking the expectation on both sides of (4) completes the proof. ∎

Lemma 5** (Bounding Variance).**

Assume that the indices in the mini-batch $\mathcal{B}_{k}$ are sampled independently from $\{1,\ldots,n\}$ with probabilities $P_{i}=\frac{L_{i}}{nL_{\mathrm{avg}}}$, then conditioned on $\bm{x}_{k}$ and $\bm{y}_{k}$, it holds that

[TABLE]

Proof.

We follow the original proofs in [31, 24] with a few modifications to take the Mahalanobis norm into account. For any $i\in\{1,\ldots,n\}$, consider the function $\phi_{i}$ defined by

[TABLE]

Then, $\phi_{i}\left(\bm{x}^{\star}\right)=\min_{\bm{x}}\phi_{i}\left(\bm{x}\right)$ since $\nabla\phi_{i}\left(\bm{x}^{\star}\right)=0$ and $\phi_{i}$ is convex. It can be checked that $\nabla\phi_{i}$ is $L_{i}$-Lipschitz, we thus have

[TABLE]

which is equivalent to

[TABLE]

Multiplying both sides of the previous inequality by $\frac{L_{\mathrm{avg}}}{L_{i}}$ and averaging from $1,\ldots,n$ gives

[TABLE]

Since $F\left(\bm{x}\right)=f\left(\bm{x}\right)+h\left(\bm{x}\right)$, by the optimality of $\bm{x}^{\star}$, $-\nabla f\left(\bm{x}^{\star}\right)\in\partial h\left(\bm{x}^{\star}\right)$. Therefore, it follows from the convexity of $h$ that $-\left\langle\nabla f\left(\bm{x}^{\star}\right),\bm{x}-\bm{x}^{\star}\right\rangle\leq h\left(\bm{x}\right)-h\left(\bm{x}^{\star}\right),$ which implies that

[TABLE]

We are now ready to bound $\operatorname{\mathds{E}}{\left\|{\Delta_{k}}\right\|_{\bm{H}^{-1}}^{2}}$, we have

[TABLE]

where the first inequality follows from the fact that $\operatorname{\mathds{E}}\left\|{X-\operatorname{\mathds{E}}X}\right\|^{2}\leq\operatorname{\mathds{E}}\left\|{X}\right\|^{2}$; we have used the inequality $\left\|{\bm{x}+\bm{y}}\right\|_{\bm{H}^{-1}}^{2}\leq 2\left\|{\bm{x}}\right\|_{\bm{H}^{-1}}^{2}+2\left\|{\bm{x}}\right\|_{\bm{H}^{-1}}^{2}$ to derive the second inequality; and the last step follows by evaluating the expectations, and using Assumption 1 and (44). This completes the proof of Lemma 5. ∎

The following definition of $\epsilon$-subgradients is very useful for analyzing how ineaxact proximal evaluations affect the convergence of the algorith.

Definition 5** ($\epsilon$-Subgradients [5]).**

Given a convex function $f:\mathbbmss{R}^{d}\to\mathbbmss{R}$ and a positive constant $\epsilon$, we say that a vector $\bm{\zeta}\in\mathbbmss{R}^{d}$ is an $\epsilon$-subgradient of $f$ at $\bm{x}\in\mathbbmss{R}^{d}$ if

[TABLE]

The set of all $\epsilon$-subgradients of $f$ at $\bm{x}$ is called the $\epsilon$-subdifferential of $f$ at $\bm{x}$, and is denoted by $\partial_{\epsilon}f\left(\bm{x}\right)$.

The following lemma characterizes the property of the $\epsilon$-differential of the function $h$, where its proof for the case of the Euclidean norm can be found in [27, Lemma 2]. We provide the proof here for completeness.

Lemma 6**.**

If $\bm{x}^{+}$ is an $\epsilon$-optimal solution to problem (40) in the sense of (5), then there exists a vector $\bm{\xi}\in\mathbbmss{R}^{d}$ such that $\left\|{\bm{\xi}}\right\|_{\bm{H}^{-1}}\leq\sqrt{2\eta\epsilon}$ and

[TABLE]

Proof.

We start by noting that $\mathbf{0}\in\partial_{\epsilon}q\left(\bm{x}^{+}\right)$ when $\bm{x}^{+}$ is an $\epsilon$-optimal solution to problem (40). Consider the function $f\left(\bm{x}^{+}\right)=\frac{1}{2\eta}\left\|{\bm{x}^{+}-\bm{t}}\right\|_{\bm{H}}^{2}$, then it can be verified that

[TABLE]

We have for convex functions $f_{1},f_{2}$ that $\partial_{\epsilon}\left(f_{1}+f_{2}\right)\left(\bm{x}\right)\subset\partial_{\epsilon}f_{1}\left(\bm{x}\right)+\partial_{\epsilon}f_{2}\left(\bm{x}\right)$ [5]. Therefore, if we let $f_{1}=f$, $f_{2}=h$, and $\bm{t}=\bm{y}-\eta\bm{H}^{-1}\bm{v}$, then $\mathbf{0}\in\partial_{\epsilon}q\left(\bm{x}^{+}\right)$. Since $\partial_{\epsilon}q\left(\bm{x}^{+}\right)\subset\partial_{\epsilon}f\left(\bm{x}^{+}\right)+\partial_{\epsilon}h\left(\bm{x}^{+}\right)$, it follows that $\mathbf{0}$ must be a sum of an element of $\partial_{\epsilon}f\left(\bm{x}^{+}\right)$ and an element of $\partial_{\epsilon}h\left(\bm{x}^{+}\right)$. Thus, there is a vector $\bm{\xi}\in\mathbbmss{R}^{d}$ such that

[TABLE]

completing the proof. ∎

Lemma 7** ([27, Lemma 1]).**

Assume that the nonnegative sequence ${u_{k}}$ satisfies the following recursion for all $k\geq 1$:

[TABLE]

where $\{S_{k}\}$ is an increasing sequence, $S_{0}\geq u_{0}^{2}$, and $\alpha_{t}\geq 0$ for all $t$. Then, for all $k\geq 1$, then

[TABLE]