Accelerated Dual Learning by Homotopic Initialization

TL;DR

This paper explores how strategic initialization in dual and non-convex optimization accelerates convergence in machine learning, supported by theoretical analysis and experimental validation.

Contribution

It introduces novel initialization strategies for dual and non-convex problems, leveraging data set properties to improve optimization speed.

Findings

Proper initialization can significantly speed up convergence.

Data correlation properties influence initial optimization progress.

Proposed heuristics are experimentally validated.

Abstract

Gradient descent and coordinate descent are well understood in terms of their asymptotic behavior, but less so in a transient regime often used for approximations in machine learning. We investigate how proper initialization can have a profound effect on finding near-optimal solutions quickly. We show that a certain property of a data set, namely the boundedness of the correlations between eigenfeatures and the response variable, can lead to faster initial progress than expected by commonplace analysis. Convex optimization problems can tacitly benefit from that, but this automatism does not apply to their dual formulation. We analyze this phenomenon and devise provably good initialization strategies for dual optimization as well as heuristics for the non-convex case, relevant for deep learning. We find our predictions and methods to be experimentally well-supported.

Tables

Table 1. Table 1: Details of the datasets.

Dataset	Size	Dimension	Boundedness ($\tau$)
a9a	32561	123	0.2
covtype	581012	54	5.5
gisette	6000	5000	100
ijcnn1	49990	22	6
w8a	49749	300	1.4
SUSY	5000000	18	1.2

Equations

$$\mathcal{Q}(\mathbf{\beta })=\frac{1}{2}{\mathbf{\beta }}^{\top }\mathbf{H}\mathbf{\beta }-{\mathbf{\beta }}^{\top }\mathbf{b},\mathbf{H}:=\frac{1}{n}{\mathbf{X}}^{\top }\mathbf{X}+\mu \mathbf{I}\in {\mathbb{R}}^{d\times d},\mathbf{b}:=\frac{1}{n}{\mathbf{X}}^{\top }\mathbf{y}\in {\mathbb{R}}^d$$

$${\nabla }_{\mathbf{\beta }}\mathcal{Q}\stackrel{!}{=}0⟺{\mathbf{\beta }}^{\ast }={\mathbf{H}}^{-1}\mathbf{b}={\left({\mathbf{X}}^{\top }\mathbf{X}+n\mu \mathbf{I}\right)}^{-1}{\mathbf{X}}^{\top }\mathbf{y}.$$

$$Q(\mathbf{\beta })$$

$${\mathbf{\beta }}^{t+1}={\mathbf{\beta }}^t-\gamma {\nabla }_{\mathbf{\beta }}\mathcal{Q}({\mathbf{\beta }}^t).$$

$${\mathbf{\beta }}^t$$

$$\mathcal{Q}(\mathbf{\beta })-{\mathcal{Q}}^{\ast }=\frac{1}{2}(\mathbf{\beta }-{\mathbf{\beta }}^{\ast })\mathbf{\Lambda }(\mathbf{\beta }-{\mathbf{\beta }}^{\ast })=\frac{1}{2}\sum _{j=1}^d{\mathbf{E}}^{\mu }[Z_j^2]({\beta }_j-{\beta }_j^{\ast }{)}^2.$$

$$\mathcal{Q}({\beta }_j^t)-{\mathcal{Q}}^{\ast }=(1-\gamma {\lambda }_j{)}^{2t}\cdot \frac{{\lambda }_j}{2}{\left({\beta }_j^0-{\beta }_j^{\ast }\right)}^2$$

$$\mathcal{Q}(0)-{\mathcal{Q}}^{\ast }=\frac{1}{2}\frac{\mathbf{E}[Y{Z}_j{]}^2}{{\mathbf{E}}^{\mu }[Z_i^2]}\le \frac{1}{2}\frac{\mathbf{E}[Y{Z}_j{]}^2}{\mathbf{E}[Z_j^2]}=\frac{1}{2}{\rho }_j^2.$$

$${\rho }_j^2=\frac{\mathbf{E}[Y{Z}_j{]}^2}{\mathbf{E}[Z_j^2]}\le \tau \mathbf{E}[Z_j^2]=\tau {\sigma }_j^2$$

$$\mathcal{Q}({\mathbf{\beta }}^t)-{\mathcal{Q}}^{\ast }\le \frac{1}{2}\left[r(\zeta )(1-\gamma \zeta {)}^{2t}+(d-r(\zeta ))\right]\tau \zeta ,r(\zeta ):=\{j:\mathbf{E}[Z_j^2]>\zeta \}$$

$${\mathcal{Q}}_{\mu }(\mathbf{\alpha })=\frac{1}{2}{\mathbf{\alpha }}^{\top }\mathbf{G}\mathbf{\alpha }-{\mathbf{b}}^{\top }\mathbf{\alpha },n\mathbf{G}=\frac{1}{\mu n}\mathbf{X}{\mathbf{X}}^{\top }+\mathbf{I},n\mathbf{b}=\mathbf{y}$$

$${\mathcal{Q}}_{\mu }(\mathbf{\alpha })$$

$${\mathcal{Q}}_{\mu }(\mathbf{\alpha })-{\mathcal{Q}}_{\mu }({\mathbf{\alpha }}_{\mu }^{\ast })=\frac{1}{2n}\sum _{i=1}^d\frac{{\lambda }_i}{\mu }{\left({\mathbf{\alpha }}_i-{\mathbf{\alpha }}_i^{\ast }\right)}^2+\frac{1}{2n}\sum _{i=d+1}^n{\left({\mathbf{\alpha }}_i-{\mathbf{\alpha }}_i^{\ast }\right)}^2$$

$${\mathbf{\alpha }}_j^{\ast \nu }={\mathbf{\alpha }}_j^{\ast \mu }(\forall j>d)\text{and}({\mathbf{\alpha }}_i^{\ast \nu }-{\mathbf{\alpha }}_i^{\ast \mu }{)}^2=n{\left(\frac{(\mu -\nu ){\sigma }_i^2}{({\sigma }_i^2+\mu )({\sigma }_i^2+\nu )}\right)}^2\mathbf{E}[Y{Z}_i{]}^2$$

$$\mathcal{Q}({\mathbf{\alpha }}^t)-{\mathcal{Q}}^{\ast }\le \frac{\tau \zeta (\nu -\mu {)}^2}{2\nu \mu }\left[r(\zeta )(1-\gamma \zeta {)}^{2t}+(d-r(\zeta ))\right],r(\zeta ):=|\{j:\mathbf{E}[Z_j^2]>\zeta \}|$$

$${\alpha }_j^{\ast }={\mathbf{u}}_j^{\top }\mathbf{y}\stackrel{\min }{⟵}\frac{1}{2n}{\alpha }_j^2-c_j{\alpha }_j,(\forall j>d)$$

$${\mathbf{\alpha }}_r^{+}={\mathbf{\alpha }}_r-{\gamma }_r{\mathcal{Q}}_r^{\prime }(\mathbf{\alpha }),$$

$${\gamma }_r^{-1}={\mathbf{G}}_{r,r}+\sum _j|{\mathbf{G}}_{r,j}|,{\gamma }_{\min }:=\underset{r}{\min }{\gamma }_r,{\gamma }_{\max }:=\underset{r}{\max }{\gamma }_r.$$

$$\mathbf{E}\Vert {\mathcal{Q}}_{\mu }^{\prime }({\mathbf{\alpha }}^{(t)}){\Vert }^2\le 2{\rho }^2\left(\frac{{\gamma }_{\max }}{{\gamma }_{\min }}\right)\Vert {\mathbf{\alpha }}_0-{\mathbf{\alpha }}^{\ast }{\Vert }^2$$

$$\mathbf{E}\left[{\mathcal{Q}}_{\mu }({\mathbf{\alpha }}^{(t)})-{\mathcal{Q}}_{\mu }^{\ast }\right]\le \frac{1}{2}{\left(1-\frac{\rho {\gamma }_{\min }}{n}\right)}^t\left({\mathcal{Q}}_{\mu }({\mathbf{\alpha }}_0)-{\mathcal{Q}}_{\mu }^{\ast }\right)+\frac{1}{2}\rho \left(\frac{{\gamma }_{\max }}{{\gamma }_{\min }}\right)\Vert {\mathbf{\alpha }}_0-{\mathbf{\alpha }}^{\ast }{\Vert }^2$$

$$\Vert {\mathbf{\alpha }}_{\nu }^{\ast }-{\mathbf{\alpha }}_{\mu }^{\ast }{\Vert }^2\le \left((\mu -\nu {)}^2/\nu \right)(d-r(\zeta ))n\tau \zeta .$$

$$\mathcal{R}(\mathbf{z})=\mathbf{E}\left[\phi ({\mathbf{x}}_i^{\top }\mathbf{z})-{\mathbf{y}}_i\left({\mathbf{x}}_i^{\top }\mathbf{z}\right)\right],$$

$${\mathbf{z}}^{+}=\mathbf{z}-{\gamma }_t{\mathcal{R}}^{\prime }(\mathbf{z})-{\eta }_t\mathbf{E}\left[y\mathbf{x}\right]$$

$$\mathcal{R}({\mathbf{z}}^{(t)})-{\mathcal{R}}^{\ast }\le c_{{\mathbf{z}}^{\ast }}\tau \varphi \left((1-\zeta /L{)}^{2t}\left(1-r(\zeta )\right)+r(\zeta )\zeta \right),$$

$$r(\zeta ):=|\{j:\mathbf{E}[Z_j^2]>\zeta \}|,c_{{\mathbf{z}}^{\ast }}:=({\mathbf{E}}_{\mathbf{x}}\left[{\phi }^{(2)}({\mathbf{x}}^{\top }{\mathbf{z}}^{\ast })\right]{)}^{-1}.$$

$${\gamma }_r={\left({\mathbf{G}}_{r,r}+\sum _j|{\mathbf{G}}_{r,j}|\right)}^{-1},{\gamma }_{\min }=\underset{r}{\min }{\gamma }_r,{\gamma }_{\max }=\underset{r}{\max }{\gamma }_r$$

$$\mathbf{E}\Vert {\mathcal{Q}}_{\mu }^{\prime }({\mathbf{\alpha }}^{(t)}){\Vert }^2\le 2{\rho }^2\left(\frac{{\gamma }_{\max }}{{\gamma }_{\min }}\right)\Vert {\mathbf{\alpha }}_0-{\mathbf{\alpha }}^{\ast }{\Vert }^2$$

$$\mathbf{E}\left[{\mathcal{Q}}_{\mu }({\mathbf{\alpha }}^{(t)})-{\mathcal{Q}}_{\mu }^{\ast }\right]\le \frac{1}{2}{\left(1-\frac{{\gamma }_{\min }\rho }{n}\right)}^t\left({\mathcal{Q}}_{\mu }({\mathbf{\alpha }}_0)-{\mathcal{Q}}_{\mu }^{\ast }\right)+\frac{1}{2}\rho \left(\frac{{\gamma }_{\max }}{{\gamma }_{\min }}\right)\Vert {\mathbf{\alpha }}_0-{\mathbf{\alpha }}^{\ast }{\Vert }^2$$

$${\mathbf{\alpha }}^{+}$$

$$=\mathbf{\alpha }-{\gamma }_r{\Psi }_r\left({\mathcal{Q}}^{\prime }(\mathbf{\alpha })-{\mathcal{Q}}^{\prime }({\mathbf{\alpha }}^{\ast })\right)$$

$$=(\mathbf{I}-{\gamma }_r{\Psi }_r\mathbf{G})(\mathbf{\alpha }-{\mathbf{\alpha }}^{\ast })+{\mathbf{\alpha }}^{\ast }.$$

$${\mathcal{L}}_{\rho }$$

Full text

1 Introduction

The standard approach to supervised machine learning is to cast it as an optimization problem with a suitable loss function and a regularizer. Learning then amounts to minimizing the regularized training risk over a chosen parametric model family. However, this view obstructs the fact that the ultimate goal is in minimizing the expected risk on unseen data and that the regularized empirical risk serves merely as a proxy for the former. Optimization for machine learning is thus by design a quest for approximate solutions. This offers learning-specific tradeoffs between statistical accuracy and computational complexity, such as early stopping [22], convex relaxation [4], data sub-sampling [7], or accepting locally optimal solutions as in deep learning [6, 21]. Whereas optimization algorithms are often analyzed in terms of their convergence to the optimum, in learning the interest should primarily be on how quickly one can find suboptimal solutions of sufficient quality (relative to the size of the data sample).

One key question in this context is what initialization strategy to use for the weights of a model, as the initial parameter choice has a huge impact on the transient phase of iterative learning algorithms. There are many examples that clearly demonstrate that initialization matters, for instance in deep learning via weight transfer [23], data-dependent initialization [15], or in avoidance of saddle points [16]. The same is true for unsupervised learning problems like matrix completion [14] or $K$-means [2]. Notably, for convex learning, the role of initialization has been somewhat neglected, mainly because of the guaranteed convergence to the global optimum. However, for massive data sets, initialization can have a huge practical impact on the scalability of an algorithm. Moreover, the (easier) convex setting also allows for a more rigorous analysis of the effect of initialization and the reasons for the slow-down caused by poor starting points. In this vein, the current paper provides a detailed analysis of convex learning, specifically of ridge regression and generalized linear models, that suggests to pre-train models with artificially increased regularization and to use this as an initialization in the spirit of homotopy or continuation methods [1, 18]. The focus of our work is on the dual problem, because it offers more flexibility in the data representation (i.e. through the use of kernels), allows for fast algorithms like stochastic coordinate descent [20] that exhibit linear convergence and is also amenable to data sharding and communication-efficient distributed implementations [13]. In particular, our method offers an yet unexplored direction to get speed-ups for sequential or parallel dual algorithms.

2 Ridge Regression

2.1 Primal Formulation

For concreteness, we perform an in-depth analysis of ridge regression [12]. Given a training set of $n$ observations $(\mathbf{x}_{i},y_{i})$, with inputs $\mathbf{x}_{i}\in\mathbb{R}^{d}$ and reponses $y_{i}\in\mathbb{R}$, we denote by $\mathbf{X}\in\mathbb{R}^{n\times d}$ the data matrix and by $\mathbf{y}\in\mathbb{R}^{n}$ the response vector. Let us assume w.l.o.g that $\mathbf{y}$ is mean zero and unit variance and that the data is centered.The ridge regression objective with regularization $\mu>0$, can be expressed as

[TABLE]

The optimal solution is explicitly given by the normal equations

[TABLE]

Note that we can interpret key quantities as expectations under the training distribution, which will subsequently become relevant. In particular note that $b_{j}=\mathbf{E}[X_{j}Y]$ is the covariance between the $j$-th feature dimension and the response. Obviously, the elements of $\frac{1}{n}\mathbf{X}^{\top}\mathbf{X}$ can be expressed as $\mathbf{E}[X_{i}X_{j}]$, i.e. they encode the empirical variances and covariances of the features. We also denote by $\mathbf{E}^{\mu}[X_{i}^{2}]=\mathbf{E}[X^{2}_{i}]+\mu$ the biased varince estimated that we get through regularization.

For the purpose of analysis and to ease the exposition, we orthogonally transform the data into an eigenfeature representation. To that extend we use the scaled SVD of $\mathbf{X}=\sqrt{n}\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}$, such that $\mathbf{X}^{\top}\mathbf{X}=n\mathbf{V}\mathbf{\Sigma}^{\top}\mathbf{\Sigma}\mathbf{V}^{\top}$ and define $\mathbf{Z}:=\mathbf{X}\mathbf{V}=\sqrt{n}\mathbf{U}\mathbf{\Sigma}$. Correspondingly the parameters are transformed via ${\bm{\beta}}\leftarrow\mathbf{V}^{\top}{\bm{\beta}}$, we arrive at the diagonalized objective

[TABLE]

2.2 Gradient Descent Analysis

Gradient descent (GD) optimizes an objective through iterative gradient updates with step size $\gamma>0$

[TABLE]

For the diagonalized quadratic objective, the iterate sequence is explicitly given as follows.

Lemma 1**.**

GD initialized at ${\bm{\beta}}^{0}$ yields the iterate sequence

[TABLE]

Corollary 2**.**

*For any $\gamma<2\lambda_{1}^{-1}$ and initial ${\bm{\beta}}^{0}$, the iterate sequence generated by GD is guaranteed to converge to ${\bm{\beta}}^{*}$ at a linear rate. *

Corollary 3**.**

Define the condition number as $\kappa:=\lambda_{1}/\lambda_{d}$. The rate of parameter convergence of $\beta_{d}$ is lower bounded by $1-\gamma\lambda_{d}>1-2/\kappa$.

Our interest is in the convergence of the objective value. We can easily (and exactly) relate distance in parameter space to suboptimality through $\mathbf{\Lambda}$.

Lemma 4**.**

Let ${\bm{\beta}}^{*}$ be the minimizer of (1) and ${\cal Q}^{*}\triangleq Q({\bm{\beta}}^{*})$. Then for any ${\bm{\beta}}$

[TABLE]

Let us compare the worst case rate we would expect based on $\kappa$ in Corollary 3 with the empirical suboptimality of the GD iterate sequence on some (randomly selected) data sets. The plots in Figure 1 show that initially (and for the relevant transient phase) the observed reduction of suboptimality is typically much better than what may be expected based on the convergence rates in parameter space. This is a striking behavior that our work aims to explain and to better exploit.

2.3 Analysis of Initialization

We study a single coordinate in the diagonalized problem.

[TABLE]

There are two factors here, one that decreases exponentially with $t$ and a constant that depends on the initialization. In an asymptotic setting, only the first would matter, leading to the well-known slow-down in convergence of GD for ill-conditioned problems, in directions where $\gamma\lambda_{j}$ becomes very small. But what if we can tolerate some suboptimality $\epsilon$? Then we could try to find an initialization such that the second term would be less than $\epsilon$. There could be two cases: (i) we devise a smart strategy to find a good ${\bm{\beta}}^{0}$ that is sufficiently close to ${\bm{\beta}}^{*}$, (ii) we set ${\bm{\beta}}^{0}=0$, but can reasonably assume that ${\bm{\beta}}^{*}$ is small. We come back to (i) and first investigate (ii) by noting that

[TABLE]

Here $\rho_{j}$ is the correlation coefficient between $Y$ and the eigenfeature $Z_{j}$. We now formulate the following hypothesis: the observed empirical sub-optimality of GD iterates is due to the fact that eigenfeatures $Z_{j}$ with small variance also have a small correlation with the response $Y$. This seems to be related to a similar assumption made in early stopping [22] and in general when using norm-based regularization: we do not want to trust features that need to be amplified a lot (low variance, but high output covariance). Hence, a reasonable feature representation should avoid encoding relevant information in such a manner. What we show here though is that such considerations not only avoid overfitting, but also accelerate GD training.

In order to subject this idea to a more formal treatment, let us introduce a suitable regularity concept. A dataset exhibits $\tau$-bounded response correlation, if for any eigen-feature $Z_{j}$ with non-zero variance

[TABLE]

With this definition, we can immediately see that we can upper bound the suboptimality of GD iterates by making use of the regularity parameter for all eigenfeatures with $\mathbf{E}[Z_{j}^{2}]\leq\zeta$, where $\zeta$ is an arbitrary threshold. Effectively, $\tau$-boundedness allows us to pay a constant approximation cost for eigenfeatures with small variance and to trade that off with the faster rates obtained for the remaining eigenfeatures. We capture the gist of this effect in the following lemma.

Lemma 5**.**

Assume that the data set has $\tau$-bounded response correlation. Then for any $\zeta\geq 0$, the suboptimality of the GD iterate sequence can be upper bounded as

[TABLE]

Proof.

Follows directly from the boundedness assumptions. ∎

Figure 2 shows that our notion of regularity seems to agree well with the empirical observations that we pointed out before. Note that regularization increases the variance of eigenfeatures and thus leads to more favorable $\tau$-bounds. However, for clarity, we wanted to capture this notion for the worst case of $\mu=0$.

3 Dual Ridge Regression

We will now consider the dual problem for ridge regression. It turns out that the initialization ${\bm{\beta}}^{0}=0$, which we analyzed for the primal in the previous section, is not advantageous for the dual. However, we can use the insight about data $\tau$-boundedness to derive an effective initialization for the dual problem.

3.1 Dual Objective

It is straightforward to derive dual objectives for convex optimization problems using Fenchel duality, e.g. see [8]. For ridge regression, we obtain the following dual objective with dual parameters ${\bm{\alpha}}$

[TABLE]

In the dual space, where typically $n\gg d$, there is an $n-d$ dimensional subspace $\text{ker}(\mathbf{X}^{\top})$, such that for the corresponding dual eigen-vectors $\mathbf{v}$, $\mathbf{v}^{\top}{\mathbf{G}}\mathbf{v}=\frac{1}{n}$. This means that we have at least $(n-d)$ orthogonal directions with a variance of $1/n$ accounting for a total variance of $(n-d)/n$. While in the primal case, there are at most $d$ terms to bound and typically $\mu\propto 1/n$, in the dual we may get a suboptimality contribution that is $\mathbf{O}(1)$, even under a $\tau$-boundedness assumption.

Similar to the primal objective, we use (to carry out the analysis, not the optimization!) the following change of variables ${\bm{\alpha}}\leftarrow\mathbf{U}^{\top}{\bm{\alpha}}$. This change of variables leads to the following diagonalization of the dual objective:

[TABLE]

and suboptimality can be written as

[TABLE]

Indeed, suboptimality can be decomposed to two terms: the suboptimality in image of the data matrix (the first sum), and the suboptimality in the kernel of the data matrix (the second sum).

Lemma 6**.**

Suppose that ${\bm{\alpha}}^{*\nu}$ is the minimizer of ${\cal Q}_{\nu}$, then for all $\nu,\mu>0$:

[TABLE]

Proof.

Follows directly from the closed form solution of minimizers. ∎

The above lemma shows that $(n-d)$ coefficients of two dual solutions with respect to different regularizer are exactly the same and the squared difference between other $d$ coefficients is propotional to $(\mu-\nu)^{2}$. Computing the minimizer ${\bm{\alpha}}^{*\nu}$ is relatively cheaper than ${\bm{\alpha}}^{*\mu}$ for $\nu\ll\mu$ as rates often depend on the strong convexity parameter. Hence we suggest initializing gradient descent for ${\cal Q}_{\mu}$ with ${\bm{\alpha}}^{*\nu}$ for a suitably chosen $\nu$.

Lemma 7**.**

Suppose that the data set has $\tau$-bounded response correlations. Then for any $\zeta\geq\frac{1}{n}$, the suboptimality of the GD iterates starting from ${\bm{\alpha}}^{*\nu}$ with step size $(n\mu)\gamma$ can be upper bounded as

[TABLE]

Proof.

Plugging the result of lemma 6 into suboptimality of the dual objective in Eq. (13) concludes the proof. ∎

The above result closely mirrors Lemma 5 with the addition of a factor that depends on the chosen $\nu$ (in relation to $\mu$). A smaller homotopic parameter $\nu$ – closer to $\mu$ – enjoys the convergence with a smaller scaling factor $(\nu-\mu)^{2}/(\nu\mu)$. However, computing the initial vector ${\bm{\beta}}_{\nu}^{*}$ is more expensive for a small $\nu$. Setting that aside, the significance of the homotopic initialization is the fact that the $(n-d)$ directions of eigenvalue $n^{-1}$ do no longer have en effect on the suboptimality. Note that in the eigensystem of Eq. (12), we can easily solve (irrespective of $\mu$)

[TABLE]

However, this is not a practical computation as the basis vectors $\mathbf{u}_{j}$ of the kernel of $\mathbf{X}$ are not known and the eigen-parameterization is not the one accessible to the algorithm. The trick of the homotopic initialization is that effectively we set $\alpha_{j}^{*}$, without having to perform the diagonalization of $\mathbf{X}^{\top}\mathbf{X}$.

3.2 Accelerated RCDM by Homotopic Initialization

In the large scale setting (when the sample size $n$ is large), the gradient step is computationally expensive. Random coordinate descent method (RCDM) [19, 20] is computationally more attractive than GD and offers a competitive convergence. A RCD-step is obtained by a coordinate-wise approximation of the gradient. More precisely, it randomly picks a random coordinate $r\in\{1,\ldots,n\}$ and updates it using the corresponding coordinate of the gradient (denoted by ${\cal Q}^{\prime}_{r}(\cdot)$) as

[TABLE]

where $\gamma_{r}$ is the coordinate-wise step size of RCDM. We ask if homotopic initialization also accelerates the convergence – up to a suboptimal solution. The convergence dependency on initialization is more subtle in a stochastic setting, where each optimization step is perturbed by the noise of the gradient approximation. Here, we theoretically prove that homotopic initialization accelerates RCDM. To this end, we provide a different convergence analysis for RCDM on the dual objective of ridge regression. Later, we will show how this result provides a better convergence rate in the objective using homotopic initialization.

Theorem 8**.**

Let ${\mathbf{G}}_{i,j}$ denote the $(i,j)$-th element of the Hessian matrix ${\mathbf{G}}$. Let the parameter vector ${\bm{\alpha}}^{(t)}$ be obtained by $t$ RCDM-steps on ${\cal Q}$, starting from ${\bm{\alpha}}_{0}$ with coordinate-wise step sizes

[TABLE]

For the above parameters, either the norm of gradient is bounded as

[TABLE]

or suboptimality is bounded as

[TABLE]

for every $\frac{1}{n}\leq\rho\leq\|{\mathbf{G}}\|$ (expectation is over the random choice of coordinates).

Suppose that $\sum_{j}|\mathbf{x}_{i}^{\top}\mathbf{x}_{j}|\leq B$ for all $i$. Using coordinate-wise step sizes of Eq. (18), the classical analysis of RCDM [19] suggests the rate $1-(n+2B/\mu)^{-1}$, while our analysis improves the rate by a factor of $\rho/n>1$ to $1-(\rho/n)(n+2B/\mu)^{-1}$. This improvement is not only with respect to the strong-convexity factor $\mu$ (e.g., as in the catalyst method [17]), but also with respect to the sample size $n$.

The acceleration of Theorem 8 is up to a suboptimal solution: either the norm of the gradient is bounded by $\mathcal{O}(\rho^{2}\|{\bm{\alpha}}_{0}-{\bm{\alpha}}^{*}\|^{2})$ or the suboptimality is at most $\mathcal{O}(\rho\|{\bm{\alpha}}_{0}-{\bm{\alpha}}^{*}\|^{2})$. Both the suboptimality bounds highly depend on the initial Euclidean distance to the minimizer. In fact, one can take advantage of the accelerated rate only if the initial distance $\|{\bm{\alpha}}_{0}-{\bm{\alpha}}^{*}\|^{2}$ is small. We observed that this distance is significantly large for initialization with the all-zero vector. Nonetheless, we theoretically bound this distance for homotopic initialization.

Lemma 9**.**

Assume that the data set has $\tau$-bounded response correlation. Then for any $\zeta$ the distance between two minimizers of the dual objective with different regularizers is bounded as

[TABLE]

Lemma 9 implies that for a smaller $\nu$, closer to $\mu$, the initial distance is smaller and hence the acceleration of homotopic initialization is up to a better suboptimality. However, convergence to ${\bm{\alpha}}^{*}_{\nu}$ is slower as it directly relates to the $\nu$. In our experiments, we observed that setting $\nu=\frac{1}{4}\sqrt{\mu}$ – which is considerably cheaper than computing the minimizer ${\bm{\alpha}}_{\mu}^{*}$– provides a significant acceleration. The coordinate-wise step sizes of Eq. (18) are quite pessimistic. In our experiments, we used larger step sizes $\gamma_{r}={\mathbf{G}}_{r,r}^{-1}$, which is equal to the coordinate-wise Lipschitz constants [19] of the dual ridge objective.

4 Generalized Linear Model

The accelerated convergence for an approximate solution also applies to generalised linear models. Given a convex differentiable smooth function $\varphi$, a generalized linear model (GLM) aims at minimizing the non-quadratic objective [9]

[TABLE]

where the expectation is over the population distribution of the data. Let $\varphi^{(k)}(\cdot)$ denote the $k$-th derivative of $\varphi(\cdot)$. The above formulation obtains logistic regression, with ${\varphi(a)=\log(1+\exp(a))}$, and ridge regression, with $\varphi(a)=a^{2}$. To take advantage of the boundedness property of a dataset, we modify the gradient descent step as

[TABLE]

where $\gamma_{t}$ and $\eta_{t}$ are two step-sizes of the modified gradient step. Indeed, we suggest to use a biased gradient step to accelerate the initial convergence. If the input vectors $\mathbf{x}_{i}$ are drawn i.i.d from a gaussian distribution, then we can obtain an accelerated convergence of the modified gradient descent on the generalized linear model. This is made precise in the next lemma.

Lemma 10**.**

Suppose that inputs are drawn i.i.d from a multivariate normal distribution with mean [math] and covariance matrix $\Sigma$, i.e. $\mathbf{x}_{i}\sim N(\vec{0},\Sigma)$ where $\|\Sigma\|=L$. Assume that $|\varphi^{(2)}(a)|\leq\phi$. If the dataset is $\tau$-bounded, then there is a step size schedule for modified gradient descent such that iterates of GD starting from zero obtains the following suboptimality bound for all $0<\zeta<L$:

[TABLE]

Although the data is assumed to be generated from a normal distribution in the last lemma, we believe that this result can be extended to an arbitrary distribution using zero-bias transformations (see [9] for more details). Furthermore, we believe that homotopic initialization can also obtain a similar acceleration, up to a suboptimal solution, on the dual objective of GLM.

5 Initialization for Deep Neural Networks

So far we have seen that if a dataset has the boundedness property, then initialization with the all-zero vector will accelerate optimization, on the primal objective, up to a sub-optimal solution. We proved this for ridge regression as well as for generalized linear models in a simplified setting. Our experiments show that non-linear features of a trained neural network provide such a representation of the data (see figure 4). This observation justifies the effectiveness of one of the most common initialization schemes of deep neural networks: Use resulting weights of a trained shallow network to initialize a deep neural network [10]. We attribute the gain of such initialization to the boundedness of features obtained by a smaller network. Suppose that we used a trained network with $N-1$ layer to initialize the first $N-1$ layers of a larger network (assume that the number of hidden units is the same as the number of output unites). If we freeze weights of the first $N-1$ layers, then optimization with respect to weights of the last layer is a convex programming task using a convex loss function, which can be modelled by a generalized linear model. Since features obtained by a trained network are $\tau$-bounded, we expect an initial acceleration of GD by initialising the last layer to all-zero. In practice, we do not need to freeze weights of early layers because back propagation naturally causes small changes in early layers due to the vanishing gradient phenomenon [11]. Furthermore, we do not need to increase the number of layers by one. In our experiments, we have observed that even a network with $N/2$ layers yields a good initialization for a network with $N$ layers (see figure 4). Since optimization of shallow networks is relatively cheaper than a deep neural network, this initialization is relatively cheap to compute. Our argument can also be extended to the layer-wise training for deep neural networks [3]. Indeed, the boundedness property is a statistical property of a representation that plays an important role in optimization for machine learning.

6 Experimental Results

6.1 Datasets and Protocol

In this section, we present our empirical results on real datasets, selected from LIBSVM library [5] (see Table 1 for more details). We measured the boundedness constant $\tau$ for these datasets, which is used in our analysis. The regularizer is $\mu=10^{-6}$ for all datasets except for Gisette, which has relatively less samples, and hence we used regularizer $\mu=10^{-3}$.

6.2 Initialization for Dual Objective

We ran an experiment to assess the advantage of the homotopic initialization on dual programming. Our experiment is on initialization of Random Coordinate Descent Method (RCDM) for optimizing the dual objective of ridge regression. Throughout all the experiments, we used the homotopic parameter $\nu=0.25\sqrt{\mu}$, which is computationally a favourable choice. We used coordinate-wise step sizes $\gamma_{r}={\mathbf{G}}_{r,r}^{-1}$, which are equal coordinate-wise Lipchitz constants and its a common choice for RCDM  [19]. Our sampling scheme of coordinates is random permutation in each epoch. Figure 3 shows the dual suboptimality, primal suboptimality, and test error through optimization. To compute the primal suboptimality, we mapped the dual parameter to the primal one using the mapping ${\bm{\beta}}^{(t)}=(n\mu)^{-1}\mathbf{X}{\bm{\alpha}}^{(t)}$ (used in [20]). This mapping obtains the primal minimizer given the dual minimizer, i.e. ${\bm{\beta}}^{*}=(n\mu)^{-1}\mathbf{X}{\bm{\alpha}}^{*}$ holds. For the test error, we computed the average of the squared loss on the test data. Overall, we observe that the homotopic initialization causes a worse initial primal and dual suboptimality compared to the initialization with all-zero vector. Nonetheless, the suboptimality decays quickly and reachs a better suboptimality compared to the initialization with zero. Indeed, the homotopic initial vector lies in the space of coarser eigen-features (features associated with larger eigenvalues) that accelerates optimization. Although the accelerated convergence is up to a suboptimal solution, this suboptimal solution achieves a test error that is comparable to the test error of the empirical minimizer on most of datasets. The acceleration, obtained by the homotopic initialization, is related to the boundedness factor $\tau$ of the dataset reported in table 1. For example, the homotopic initialization on gisette, which has a large $\tau$, obtains relatively less gain.

Although our analysis of homotopic initialization was limited to the ridge regression problem, our experiments on dual SVM show the same behaviour for this initialization. Our results on dual SVM are included in the appendix.

6.3 Initialization for neural networks

We ran an experiment to highlight the role of boundedness of a representation in the convergence speed of gradient descent on neural networks. In this experiment, we train a multilayer perceptron (MLP) with 10 hidden layers and 100 hidden units in each layer. We used two datasets MNIST and CIFAR with 50000, and 20000 samples, respectively. Here, we compared the boundedness of features obtained from the last layer of the network before and after training. Our observations show that a trained neural network provides boundedness (see figures 4.a and 4.b). Based on this observation, we use a layer-wise initialization strategy. We trained a MLP with 5 layers and the same number of hidden units. Then we initialized the first 5 layers of the main network (with 10 layers) by the trained weights of the smaller network and we set the rest of weights to zero.We compared the convergence of GD with step size 0.1 using these two different initialisation schemes. The layer-wise initialization, which yields a representation with the boundedness property, significantly accelerates the initial convergence of GD.

Acknowledgments

We would like to thank Aurelien Lucchi, Octavian Ganea, and D nner Celestine for helpful discussions.

7 Appendix

Here, we present the proof of lemmas and theorems of the manuscript with our experiments on the SVM problem. We start with the analysis of random coordinate descent method, then we provide the proof of the claim on generalized linear model. Finally, we represent the experimental result.

7.1 Analysis of RCDM

We proved that homotopic initialization accelerates convergence of gradient descent, up to a suboptimal solution, on the dual objective. In this section, we extend the result to the random coordinate descent method. For gradient descent, iterates could be tracked in a closed form, and the role of the initialisation was closely reflected in the convergence. The analysis of random coordinate descent, however, is more involved due to the noisy estimation of the gradient. In the next lemma, we decompose the convergence bound of RCDM to two terms: a term that depends on the initial sub-optimality on the dual objective, the second term is determined by the initial Euclidean distance to the dual minimizer. Later, we will use this result to prove that homotopic initialization provides an acceleration up to a suboptimal.

Theorem** (8).**

Let ${\mathbf{G}}_{i,j}$ denote the $(i,j)$-element of the Hessian matrix ${\mathbf{G}}$. Let the parameter vector ${\bm{\alpha}}^{(t)}$ is obtained by $t$ RCDM-steps on ${\cal Q}$ starting from ${\bm{\alpha}}_{0}$ with coordinate-wise step sizes

[TABLE]

For the obtained parameter, either the norm of gradient is bounded as

[TABLE]

or suboptimality is bounded as

[TABLE]

for every $\frac{1}{n}\leq\rho\leq\|{\mathbf{G}}\|$ (expectation is over the random choice of coordinates).

Proof.

We rewrite the RCDM update in a form that facilitates our future analysis. Consider matrix $\Psi_{r}$ whose $r$-th diagonal element is one and remaining elements are zero; using this matrix, RCDM step can be written as

[TABLE]

Note that we used the optimality condition of ${\bm{\alpha}}^{*}$ in the second step. We skipped the subscript $\mu$ for ${\cal Q}_{\mu}$ throughout the proof. To obtain the desired convergence guarantee, we decompose the Hessian matrix, of the dual objective, as ${\mathbf{G}}={\cal L}_{\rho}+{\cal S}_{\rho}$ where matrices ${\cal L}_{\rho}$ and ${\cal S}_{\rho}$ are obtained from SVD decomposition of the data matrix:

[TABLE]

The second smallest eigenvalue of the ${\cal L}_{\rho}$ is $\rho$, the spectral norm of ${\cal S}_{\rho}$ is bounded by $\rho$, and ${\cal L}_{\rho}{\cal S}_{\rho}=0$. Using the above decomposition, the expected suboptimality after one RCD-step can be written as

[TABLE]

We prove that the first term, which depends on directions with large eigenvalues, decays in a favourable rate (in lemma 11) and the second term is bounded (by a factor of the threshold $\rho$).

Lemma 11**.**

RCMD with coordinate-wise step sizes

[TABLE]

guarantees the decrement

[TABLE]

as long as

[TABLE]

Proof.

Consider the compact notations $\Delta={\bm{\alpha}}-{\bm{\alpha}}^{*}$ and $\Delta_{+}={\bm{\alpha}}^{+}-{\bm{\alpha}}^{*}$. We rewrite the expected suboptimality after one setp of RCDM as

[TABLE]

where terms $T_{1}$ and $T_{2}$ are formulated as

[TABLE]

The first term is bounded as:

[TABLE]

In steps [1]–[3], we used the fact that $\mathbf{E}\left[\Psi_{r}\right]=\mathbf{I}/n$, ${\mathbf{G}}={\cal L}_{\rho}+{\cal S}_{\rho}$, and the second smallest eigenvalue of the matrix ${\cal L}_{\rho}$ is greater than $\rho$, respectively. It remains to show that $T_{2}$ is negative:

[TABLE]

where we used following facts in each step:

• [1]:

Note that ${\cal L}_{\rho,r}$ denotes the diagonal element $r$ of matrix ${\cal L}_{\rho}$. Matrix $\Psi_{r}$ has only one non-zero element: the $r$-th diagonal which is one, hence $\Psi_{r}{\cal L}_{\rho}\Psi_{r}={\cal L}_{\rho,r}\Psi_{r}$.

• [2]:

Recall $({\cal L}={\mathbf{G}}-{\cal S})\rightarrow({\cal L}_{\rho,r}={\mathbf{G}}_{r}-{\cal S}_{\rho,r})$. Since ${\cal S}_{\rho}$ is positive definite, its diagonal elements are positive. This concludes that ${\cal L}_{\rho,r}\leq{\mathbf{G}}_{r}$.

• [3]:

The choice of step sizes leads to $\gamma_{r}{\mathbf{G}}_{r}\leq 1/2$.

• [4]:

We know that $\mathbf{E}_{r}\left[\Psi_{r}\right]=\mathbf{I}/n$ according to the definition of matrix $\Psi_{r}$.

• [5]:

Matrix ${\cal L}_{\rho}$ and ${\cal S}_{\rho}$ are orthogonal.

• [6]:

As long as the norm of the gradient is sufficiently large (the condition in Eq. (35)), this inequality holds. To prove this, we write the norm of the gradient as:

[TABLE]

Now the condition of Eq. (35) implies

[TABLE]

∎

The last lemma implies that the convergence of $\mathbf{E}\left[\|{\bm{\alpha}}^{+}-{\bm{\alpha}}^{*}\|_{{\cal L}_{\rho}}\right]$ is dominated by $\rho$ up to a suboptimal solution in term of the norm of the gradient. It remains to prove that the second term of suboptimality, i.e. $\mathbf{E}\left[\|{\bm{\alpha}}^{+}-{\bm{\alpha}}^{*}\|_{{\cal S}_{\rho}}\right]$, scales with a factor of $\rho$:

[TABLE]

Let parameter vector ${\bm{\alpha}}^{(t)}$ denote the obtained parameter vector after $t$ RCDM steps. We know that RCMD with step size of $\gamma_{r}\leq({\mathbf{G}}_{r,r})^{-1}$ monotonically convergences to the minimizer in terms of the objective value. Furthermore, $(1/n)$-strong convexity of the dual function implies that $\|{\bm{\alpha}}^{(t)}-{\bm{\alpha}}^{*}\|^{2}\leq n({\cal Q}^{(t)}-{\cal Q}^{*})$. Therefore, $\|{\bm{\alpha}}^{(t)}-{\bm{\alpha}}^{*}\|$ asymptotically convergences to zero; however, this convergence is not monotone. In other words, the distance might diverge. In the next lemma, we bound this distance by a factor of the initial distance, i.e. $\|{\bm{\alpha}}_{0}-{\bm{\alpha}}^{*}\|^{2}$.

Lemma 12**.**

The distance of RCD iterates ${\bm{\alpha}}^{(t)}$ from the minimizer are bounded by a constant factor of the initial distance from the minimizer:

[TABLE]

for the choice of step sizes in Eq. (26).

Proof.

We introduce the diagonal matrix $\Lambda$ whose $r$ diagonal element is the coordinate-wise step size $r$. The definition of $\Lambda$ implies $\gamma_{\min}\mathbf{I}\preceq\Lambda\preceq\gamma_{\max}\mathbf{I}$. Furthermore, the choice of step sizes of Eq. (26) ensures matrix $\Lambda-{\mathbf{G}}$ is diagonal dominant, therefore it is positive definite. Using these facts, we prove that $\mathbf{E}\|{\bm{\alpha}}^{(t)}-{\bm{\alpha}}^{*}\|^{2}_{\Lambda}$ does not diverge. Some straightforward algebra leads to

[TABLE]

where $T_{1}$ is a negative definite matrix:

[TABLE]

where steps [1]–[3] are obtained by the choice of step size, $\mathbf{E}\left[\Psi_{r}\right]=\mathbf{I}/n$, and $0\preceq\Lambda-{\mathbf{G}}$, respectively. Using this result, we bound the desired distance as

[TABLE]

Step [1] is obtained by recursion on $t$. Step [2] is derived by the definition of the matrix $\Lambda$. ∎

Using last two lemmas, we prove the theorem:

[TABLE]

where we used following facts:

• [1]:

This step is obtained from spectral decomposition of the hessian matrix (see Eq. (32)).

• [2]:

Recursion on the result of lemma 11 yields the inequality.

• [3]:

This inequality is a direct result of lemma 12.

∎

The suboptimality bound of the last theorem provides a fast convergence rate up to a suboptimal solution that is dominated by the initial distance $\|{\bm{\alpha}}_{0}-{\bm{\alpha}}^{*}\|^{2}$. In the next lemma, we bound this distance using the homotopic parameter $\nu$.

Lemma** (9).**

Assume that the data set is $\tau$-bounded. Then the distance between two minimizers of the dual objective with different regularizers is bounded as

[TABLE]

Proof.

Lemma 6 with some straight forward algebra conclude the result. ∎

7.2 Generalized Linear Model

The gradient of GLM

Recall the objective of GLM:

[TABLE]

Suppose that the $\mathbf{x}$ are from a multivariate normal distribution, i.e. $\mathbf{x}\sim N(\vec{0},\Sigma)$. Then one can use integration by parts (Stein’s lemma [9]) to write the gradient of the above objective as

[TABLE]

where $\varphi^{(2)}$ denote the second derivate of $\varphi(.)$. The Hessian matrix of ${\cal R}$ highly depends on the covariance matrix of the distribution of input:

[TABLE]

Using the above the Hessian and the gradient, we provide a useful expression of the gradient that facilitates the our convergence analysis.

Lemma 13**.**

There exists constants $\xi_{1}$, and $\xi_{2}$, which depend on $\mathbf{z}$, and $\mathbf{z}^{*}$ such that

[TABLE]

holds.

Proof.

According to mean-value theorem, there is a $\bar{\mathbf{z}}=(1-\theta)\mathbf{z}+\theta\mathbf{z}^{*},\theta\in\left(0,1\right)$ such that

[TABLE]

where constants are

[TABLE]

The Eq. 59 with the optimality condition ${\cal R}^{\prime}(\mathbf{z}^{*})=0$ conclude the proof for $\xi_{1}=c_{2}/(1-c_{3})$, and $\xi_{2}=c_{3}/(1-c_{3})$. ∎

A biased gradient step

We suggest our modified gradient step as

[TABLE]

where constant $\gamma_{t}$ and $\eta_{t}$ are two step sizes. In the next lemma, we prove that the convergence of such a modified gradient descent depends on the covariance matrix of the distribution.

Lemma 14**.**

Let iterate $\mathbf{z}^{(t)}$ is obtained by $t$ GD steps on GLM. There is a schedule for $\gamma_{t}$ and $\eta_{t}$ such that modified GD steps can be written as

[TABLE]

where $L=\|\Sigma\|$.

Proof.

We prove the above bound by induction on $t$.

[TABLE]

The above equation with step size $\gamma_{t}=(L\xi_{1})^{-1}$ and $\eta=-\gamma_{t}\xi_{2}$ completes the proof.

holds. Replacing the above equation in the convergence bound with the step size $\gamma_{t}=(L\mathbf{E}_{\mathbf{x}}\left[\varphi^{(2)}(\mathbf{x}^{\top}\bar{\mathbf{z}}^{(t)})\right])^{-1}$ concludes the proof. ∎

Although the last lemma does not specify step sizes, one can find this schedule by line-search. The result allows us to track iterates of GD in a closed form, which provides an accelerated convergence up to a suboptimal solution for GLM – similar to ridge regression. The next lemma proves this.

Lemma** (10).**

Suppose that inputs are drawn i.i.d from a multivariate normal distribution with mean [math] and covariance matrix $\Sigma$, i.e. $\mathbf{x}_{i}\sim N(\vec{0},\Sigma)$ where $\|\Sigma\|=L$. Assume that $|\varphi^{(2)}(a)|\leq\phi$. If the dataset is $\tau$-bounded, then there is a step size schedule for modified gradient descent such that iterates of GD starting from zero obtains the following suboptimality bound for all $0<\zeta<L$:

[TABLE]

Proof.

Using the mean-value theorem, there is $\bar{\mathbf{z}}^{(t)}=\theta\mathbf{z}^{(t)}+(1-\theta)\mathbf{z}^{*},\theta\in\left[0,1\right]$ such that

[TABLE]

It remains to relate the $\mathbf{z}^{*}$ to the boundedness assumption. To this end, we use optimality condition of $\mathbf{z}^{*}$ as

[TABLE]

The above result, which is also provided in [9], shows that the minimizer of the GLM can be obtained by scaling the minimizer of ridge regression. The convergence bound of Eq. (72) is also scaled convergence of the gradient descent on ridge regression. Plugging the minimizer into the Eq. (72) with the boundedness assumption completes the proof. ∎

7.3 Experiments