Title: TRANSFERRING OPTIMALITY ACROSS DATA DISTRI-BUTIONS VIA HOMOTOPY METHODS

Abstract: Homotopy methods, also known as continuation methods, are a powerful mathematical tool to efficiently solve various problems in numerical analysis. In this work, we propose a novel homotopy-based numerical method that can be used to gradually transfer optimized parameters of a neural network across different data distributions. This method generalizes the widely-used heuristic of pre-training parameters on one dataset and then fine-tuning them on another dataset of interest. We conduct a theoretical analysis showing that, under some assumptions, the homotopy method combined with Stochastic Gradient Descent (SGD) is guaranteed to converge in expectation to an r θ -optimal solution for a target task when started from an expected r θ -optimal solution on a source task. Empirical evaluations on a toy regression dataset and for transferring optimized parameters from MNIST to Fashion-MNIST and CIFAR-10 show substantial improvement of the numerical performance over random initialization and pre-training.

1. INTRODUCTION:

Homotopy methods (Allgower & Georg, 1980) , also known as continuation methods, are a powerful mathematical tool to efficiently solve various problems in numerical analysis (e.g., Tran-Dinh et al. (2012) , Zanelli et al. (2019) ). The core idea consists in sequentially solving a series of parametric problems, starting from an easy-to-solve problem and progressively deforming it, via a homotopy function, to the target one. Homotopy methods are suitable to solve complex non-convex optimization problems where no or only little prior knowledge regarding the localization of the solutions is available. In addition, in contrast to state-of-the-art algorithms in deep learning (e.g., Bottou (2010), Duchi et al. (2011) , Kingma & Ba (2015) ), these methods often achieve global convergence guarantees by only exploiting local structures of the problem. Concepts, such as curriculum-learning and warm-starting, that are related to different degrees to homotopy methods, have been explored both in the deep learning (e.g., Gulcehre et al. (2016) , Mobahi (2016) , Gulcehre et al. (2017) ) and in the reinforcement learning (e.g., Narvekar (2017) ) communities.

In this work, we propose a novel homotopy-based numerical method to transfer knowledge regarding the localization of a minimizer across different task distributions in deep learning. This method gradually tracks a neural network's (close-to-)optimal parameters from one data distribution to another one via the homotopy method (Allgower & Georg, 1980) and can be interpreted as a generalization of the very common heuristic of fine-tuning a pre-trained network. After discussing related work (Section 2) and background on homotopy methods (Section 3), our contributions are as follows:

1. We provide a general theoretical analysis of the homotopy method when using SGD as an iterative solver, proving that under some local assumptions it tracks in expectation an r θ -optimal solution from the source task to the target task (Section 4).

2. We introduce homotopy functions for transferring optimality across data distributions for supervised regression and classification tasks (Section 5).

3. For a toy regression dataset and for transferring optimized parameters from MNIST to Fashion-MNIST and from MNIST to CIFAR-10, we show that our method obtains up to two orders of magnitude better numerical performance than random initialization and substantial improvement of the numerical performance over pre-training (Section 6).

2. RELATED WORK:

Deep neural networks have led to establish a new state-of-the-art in many applications. Despite their great success and the many theoretical studies that have been published in the last years (e.g., Balduzzi et al. (2017) , Li et al. (2018) , Feizi et al. (2018) , Kunin et al. (2019)) , training these deep models remains a big challenge. Various stochastic optimization algorithms (e.g., Duchi et al. (2011) , Kingma & Ba (2015) , Reddi et al. (2018) ) and initialization heuristics (e.g., Daniely et al. (2016) , Klambauer et al. (2017) , Hanin & Rolnick (2018) ) have been recently suggested in order to improve and speed up the training procedure. We now briefly discuss the state-of-the-art deep learning optimization techniques and initialization strategies that are most related with the proposed homotopy-based method, drawing connections with existing and ongoing research works in the field.

Curriculum Learning. First introduced by Bengio et al. (2009) and then extended in different works (e.g., Graves et al. (2017) , Weinshall et al. (2018) , Hacohen & Weinshall (2019) ), curriculum learning can also be listed among the optimization heuristics proposed to alleviate the complexity of solving high dimensional and non-convex problems. In particular, taking inspiration from the fact that humans and animals learn "better" when exposed to progressively more complex situations in an organized manner, curriculum learning techniques guide the training by starting with "easy-tolearn" samples and progressively introducing more "complex-to-learn" ones. This guided learning process can also be rephrased in a homotopy-like fashion (see Algorithm 1) as solving a sequence of optimization problems where the target training distribution gradually changes from considering only the "easy" examples to the full original training distribution.

Meta-Learning and Transfer-Learning. Due to the massive amount of computational resources required by the development of modern deep learning applications, the community has started to explore the possibility of re-using learned parameters across different tasks, leading to the development of many new transfer-learning (e.g., Rohrbach et al. (2013) , Wang & Schneider (2014) , Cui et al. (2019) ) and meta-learning (e.g., Schmidhuber (1987 ), Hochreiter et al. (2001 ), Finn et al. (2017) , Zintgraf et al. (2019) ) algorithms. The simplest way to transfer knowledge across different tasks consists in using warm-start initialization. This heuristic is amply used in computer vision applications, where it is also known as the fine-tuning technique (e.g., Krizhevsky et al. (2012) , Yosinski et al. (2014) , Reyes et al. (2015) , Käding et al. (2016) ). So far, there is no rigorous explanation of why and when fine-tuning works. However, numerous empirical evaluations on different benchmarks show that warm-starting the parameters of deep models often leads to faster convergence and better generalization than using random initialization.

3. BACKGROUND:

In this work, we will focus on solving problems of the form

EQUATION

where J : R d → R is our target objective function and θ * is a minimizer. Problems as described in (1) arise, for instance, in classification and regression scenarios.

In the following section we briefly review the main concepts of homotopy and continuation methods, which the proposed technique to solve problem (1) is based on.

3.1. HOMOTOPIC FUNCTIONS AND CONTINUATION METHODS FOR OPTIMIZATION:

Given two topological spaces Z and Y , a homotopy is a continuous deformation between two continuous functions g, f : Z → Y that fulfills certain properties. We can formalize this concept with the following definition Definition 3.1. Let g, f : Z → Y be continuous maps on the topological spaces Z, Y . A homotopy from g to f is a continuous function

EQUATION

If such function H exists, g is said to be homotopic of f , and this relation is denoted by g f .

It is straightforward to show that, A ⊆ R n being a convex set, any two continuous maps g, f : Z → A are homotopic (see (Suciu, 2016) for a derivation). From this fact it follows that any two continuous and real functions are homotopic. See Figures 4a-4b in the appendix for a graphical representation of two different homotopy maps between the probability density functions of two Gaussian distributions, where λ ∈ [0, 1] denotes the homotopy parameter. See also Section A in the appendix for details on some of the main properties of homotopic functions.

Continuation methods (also known as homotopy methods) are a widely used mathematical tool to solve complex non-convex optimization problems where no or only very limited prior knowledge regarding the localization of optimal solutions is available (see (Allgower & Georg, 1980) for a full characterization of continuation methods). The core idea of a homotopy approach consists in defining a homotopy function H(θ, λ) with λ ∈ [0, 1] such that H(θ, 0) = J 0 (θ) is a trivial to optimize smooth map (or a smooth map of which a surrogate θ 0 of an optimal solution is available) and H(θ, 1) = J(θ) is our target objective function. Instead of directly addressing problem (1), we approximately and sequentially solve γ > 0 parametric optimization problems of the form

EQUATION

for increasing values of the parameter λ i for i = 1, . . . , γ and warm-starting each problem with the previously derived approximate solution. Conceptually, Algorithm 1 describes the basic steps of a general homotopy algorithm. Under appropriate assumptions, if the increment ∆λ is sufficiently small, then the iterative procedure in Algorithm 1 will converge to a neighborhood of an optimal solution of the target objective J that depends in some sense on the number of iterations k > 0 performed (Allgower & Georg, 1980) . Many different variations of Algorithm 1 exist. In particular,

Algorithm 1 A Conceptual Homotopy Algorithm 1: θ 0 ≈ θ * 0 ∈ arg min θ H(θ, 0) 2: γ > 0 , γ ∈ Z 3: λ 0 = 0, ∆λ = 1/γ 4: k > 0 , k ∈ Z 5: for i = 1, . . . , γ do 6: λ i ← λ i-1 + ∆λ 7: procedure θ i ←ITERATIVESOLVER(θ i-1 , k, H(θ, λ i )) 8: return θ γ

different update schemes for the homotopy parameter can be adopted (e.g., geometric or sublinear rate of increase), various iterative solvers can be used under distinct and specific assumptions, and, finally, also diverse levels of approximation for the solutions θ * i can be considered, i.e. different k values.

Before going into the details of two concrete formulations of the conceptual homotopy method outlined in Algorithm 1 (see Section 5) when applied to transfer optimality knowledge in regression and classification scenarios, we provide a general theoretical analysis in a simplified setting.

4. THEORETICAL ANALYSIS:

In this section, we provide a local theoretical analysis of homotopy methods when Stochastic Gradient Descent (SGD) (Bottou, 2010) is used as iterative solver in Algorithm 1. The locality of the analysis consists in the definition of hyperspheres of radius B ≥ 0 around the optimal solutions of each homotopy problem H(θ, λ i ) where it is possible to exploit certain structures of the problem. In particular, we approximately and sequentially solve γ > 0 unconstrained optimization problems of the form θ * i ∈ arg min

EQUATION

where H(θ, λ i ) fulfills the assumptions described in Section 4.1 and λ i ∈ [0, 1]. Let θ i be an approximate solution of the problem associated with parameter λ i derived by applying k > 0 iterations of SGD (in the limit, k = 1) and also the starting point for the problem associated with parameter λ i+1 , ∀i = 1, . . . , γ -1. In addition, let θ 0 denote an approximate solution for the source task, i.e. λ 0 = 0, that is used as initial point for the problem associated with λ 1 . In this section we characterize the maximum allowed variation of the homotopy parameter in order for the method to able to track in expectation an r θ -optimal solution from source to target task.

4.1. ASSUMPTIONS:

We now expose the fundamental assumptions for our general local theoretical analysis on which all the derivations in Sections 4.2 and 4.3 rely. In addition, throughout the analysis the -functions in (3) are implicitly assumed to be differentiable in θ. We start by giving the definition of the regions around the optimal solutions of the homotopy problems where the analysis is conducted. Definition 4.1. Given θ * i and B ≥ 0, let B B,θ * i be the following set of vectors B B,θ * i := {θ s.t. θ -θ * i ≤ B} , ∀i = 0, . . . , γ . Assumption 4.2 (local L-smoothness). Assume that there exists a constant L > 0 such that

EQUATION

Corollary 4.2.1. If H is locally L-smooth in θ, then the following inequality holds

EQUATION

Proof. See Lemma 1.1 in (Gower, 2018) for a proof.

Assumption 4.3 (local µ-strong convexity). Assume that there exists µ > 0 such that

EQUATION

Assumption 4.5 (local bounded "variance"). Let g( θ, λ i ) denote an unbiased estimate of the gradient ∇ θ H( θ, λ i ). Assume that there exists a constant C ≥ 0 such that the following bound on the expected squared norm of the estimate of the gradient holds

EQUATION

Remark 4.6. Assumption 4.5 is standard for proving error bounds on SGD iterates (see (Schmidt, 2014) ). In addition, notice that, since

E g( θ, λ i ) 2 = Var g( θ, λ i ) + E g( θ, λ i ) 2 ,

the C constant is proportional to the variance and the squared expected value of the norm of the gradient estimate. Therefore, it decreases when the iterates approach a minimizer and by reducing the noise in the estimate of the gradient. In the limit (i.e. exact gradient and convergence to a minimizer), C = 0.

Recall that θ * (λ i ) ≡ θ * i . Assumption 4.7 (strong regularity). Assume that there exists δ > 0 such that the following inequality holds θ * (λ i+1 ) -θ * (λ i ) ≤ δ|λ i+1 -λ i | , ∀i = 0, . . . , γ -1. Remark 4.8. Assumption 4.7 follows directly from the application of the Implicit Function Theorem by introducing some milder assumptions on the problem structure (see Lemma 2.1.8 in (Allgower & Georg, 1980) ).

4.2. FUNDAMENTAL THEORETICAL PRELIMINARIES:

Before proceeding with the main theoretical contributions, we extend the existing results in the literature on global error bounds for the iterates of Stochastic Gradient Descent such that they can be applied when the underlying assumptions are only required to hold locally. The derived local error bounds for SGD iterates are used in Proposition 4.11 and Theorem 4.12. Proposition 4.9. Let θ i ∈ B B,θ * i be the starting point for the problem described in (3), and let θ i := θ i,0 and θ i+1 := θ i,k denote the iterate after k > 0 SGD steps, where an SGD step is defined as

θ i,k = θ i,k-1 -αg(θ i,k-1 , λ i ) .

Under Assumptions 4.2-4.5 and by setting the batch size 0 < M ≤ N to a value such that

(N -M ) N ≤ (1-κ d ) 2αν B with κ d = (1 -αµ)

and the learning rate α to a constant value such that 0 < α ≤ min 1 2µ , 1 L , the following error bound on the iterates holds

EQUATION

Proof. See Section D in the appendix.

Remark 4.10. The expectation in (10) is taken w.r.t. all the random variables, i.e. estimates of the gradients and initial point θ 0 , involved in the optimization procedure up to the current i + 1 iteration of the algorithm.

4.3. MAIN THEORETICAL CONTRIBUTIONS:

Under the considered assumptions and by exploiting the previously derived results on local error bounds for SGD iterates, we show that, if the approximate solution θ i for the problem with parameter λ i is "sufficiently close" to a minimizer θ * i in expectation, i.e. E θ i -θ * i 2 ≤ r 2 θ , then, for a "sufficiently small" change in the homotopy parameter, the same vicinity to a minimizer θ * i+1 is preserved in expectation for the approximate solution θ i+1 of the problem with parameter λ i+1 , i.e.

E θ i+1 -θ * i+1 2 ≤ r 2 θ .

In particular, with Theorem 4.12 we characterize the maximum allowed variation of the homotopy parameter based on the properties of the parametric problems and the convergence characteristics of the adopted iterative solver, i.e. rate of convergence and number of iterations.

First, in order to apply the results derived in Theorem 4.12, given a realization of θ i ∈ B B,θ * i , we have to derive the conditions on

θ i -θ * i such that θ i -θ * i+1 ≤ B.

In addition, we derive the necessary conditions in order to apply these results recursively across the iterations of Algorithm 1. Proposition 4.11.

Let θ i ∈ B B,θ * i and |λ i -λ i+1 | ≤ , with 0 ≤ ≤ B δ . If θ i -θ * i ≤ B -δ , then θ i -θ * i+1 ≤ B. Moreover, let κ d = (1 -αµ) and assume that (N -M ) N ≤ (1 -κ k d )(1 -κ d )B 2αν ,

and

≤ 1 δ (1 -κ k d )B - (N -M ) N 2αν (1 -κ d )

.

Then, after applying k iterations of SGD, we obtain that

θ i+1 -θ * i+1 ≤ B -δ .

Proof. See Section E.1 in the appendix.

See Figure 9 in the appendix for a graphical representation of the results derived in Proposition 4.11, where the continuous and dashed lines are used to represent the circles of radius B and B -δ , respectively. Theorem 4.12. Consider Algorithm 1 with Stochastic Gradient Descent as solver and let k > 0 be the number of iterations, 0 < α ≤ min 1 2µ , 1 L be the step size and 0 < M ≤ N be the batch size such that

(N -M ) N ≤ (1 -κ k d )(1 -κ d )B 2αν ,

where

EQUATION

EQUATION

)

the following inequality holds

EQUATION

Proof. See Section E.2 in the appendix.

The results derived in Theorem 4.12 show that the homotopy method used in combination with SGD allows to track in expectation an r θ -optimal solution across the parametric problems for "small enough" variations of the homotopy parameter, i.e. ∆λ ≤ ˜ . Notice that r θ can potentially be smaller than B -δ and has to be bigger than the radius of the "noise-dominant" hypersphere centered at the minimizers, i.e. r 2 θ ≥ αC 2 2µ . In particular, by exploiting the local structure of the parametric problems we derive the maximum allowed variation of the homotopy parameter across the iterations of Algorithm 1. The derived upper bound is inversely proportional to the strong regularity constant δ and depends on the number of iterations k performed with SGD, such that the more iterations we perform on each parametric problem the more we are allowed to change the homotopy parameter. Finally, notice that these results can be applied recursively across the parametric problems.

5. TRANSFERRING OPTIMALITY VIA HOMOTOPY METHODS:

In this section we describe a possible application of homotopy methods to solve supervised regression and classification tasks. We address the case where deep neural networks are used as models. We start by introducing the problem framework of supervised learning and then we propose two different homotopy functions for the regression and classification scenarios, respectively.

5.1. PROBLEM FORMULATION:

Despite the generality of the proposed methodology, in this work we specifically address the supervised learning framework, and, in particular, when the predictive model is constituted by a deep neural network f (x; θ) parameterized by θ ∈ R d .

In the supervised learning scenario, independently from the type of task t, we typically dispose of a training set D t consisting of N pairs of examples (x j , y j ). The goal of the learning process is to find a value of θ that minimizes an objective function which measures the discrepancy between the outputs produced by the network ŷ = f (x; θ) and the target outputs y. In particular, the learning process consists in minimizing the following empirical objective function

EQUATION

whose non-convexity originates from the high non-convexity of our model f .

In the classical setting, J is chosen based on the KL divergence between the target data distribution Q x,y , with density q x,y = q(y|x)q(x), and the learned data distribution P x,y (θ), with density p x,y = p(y|x; θ)q(x), where p(y|x; θ) is modeled via a neural network, (Goodfellow et al., 2016) . With the appropriate approximations, this leads to the following form for the objective function

EQUATION

)

5.2 HOMOTOPY FUNCTIONS ACROSS DATA DISTRIBUTIONS Finding a value of θ that attains a local minimum of the objective function in ( 14) is often a hard optimization task, given the high dimensionality and non-convexity of the problem. In addition, prior knowledge regarding the localization of the solutions is rarely available. The complexity of minimizing such functions also depends in some non-trivial way on the task distribution Q x,y that is addressed (e.g., Ionescu et al. ( 2016), Zendel et al. (2017) ). For some tasks, convergence to a good approximate solution is achieved after a few epochs, while for other tasks, orders of magnitude more iterations are required to reach the neighborhood of a solution. In this perspective, different heuristics have been recently proposed in the attempt of re-using across different data distributions the prior knowledge gained from approximately solving the learning problem associated with a certain task. The question whether we could exploit easy-to-solve or already-solved tasks to speed up and improve the learning of unsolved hard tasks arises. The method we propose in this paper addresses this question and attempts to do so by using a rigorous and well-established mathematical framework, with the goal of speeding up the learning process in presence of hard-to-solve tasks.

In the perspective of homotopy methods, this goal can be achieved under some assumptions by defining a homotopy transformation between starting and target tasks and by following the procedure described in Algorithm 1. Despite the flexibility and generality of the method, with this work we only focus on homotopy deformations across different task distributions, but similar transformations can be applied in numerous different manners that are also worth exploring, e.g., progressively modifying the architecture of the network or the weights of the objective function terms.

Let s be the source task with training data D s of pairs (x s , y s ) ∼ Q xs,ys whose good approximate solution θ * s for the minimization of the objective in ( 14) is available (or cheaply computable), and let t denote the target task with training data D t of pairs (x t , y t ) ∼ Q xt,yt whose conditional distribution we aim to learn. We propose two different homotopy deformations from task s to task t for regression and classification, respectively.

5.2.1. SUPERVISED REGRESSION:

In the supervised regression scenario, by modeling the density of the conditional learned distribution as p(y|x; θ) = N y; f (x, θ), σ 2 I and using the approximate KL divergence objective function described in (15), we recover the mean squared error as minimization criterion. The proposed homotopy deformation is based on the following equations

EQUATION

EQUATION

)

Notice that the transformation described in ( 16) preserves the unimodality of the conditional distribution (see caption of Figures 4a and 4b in the appendix), and, when used in combination with the objective function defined in Equation ( 15), leads to the minimization w.r.t. θ of

EQUATION

)

See Figure 6a in the appendix for a graphical representation of this homotopy deformation when applied to gradually transform a one-dimensional sine wave function with a frequency of 1 radian into a one-dimensional sine wave function with a frequency of 137 radians. A downside of this homotopy deformation is that the same support for x is required (the absence of the subscripts s and t on x stands to indicate that the same realization for x s and x t has to be considered). Alternatively, it is possible to approximate (16) by using a Gaussian filter (see Figure 6b and Section B in the appendix).

5.2.2. SUPERVISED CLASSIFICATION:

In the case of supervised classification, by modeling the density of the conditional learned distribution as p(y|x; θ) = M ultinoulli(y; f (x; θ)), and using the approximate KL divergence objective function described in (15), we recover the cross-entropy loss function, (Goodfellow et al., 2016) . A possible homotopy deformation for the classification case consists in applying the following transformations

EQUATION

y λ |x λ = (1 -λ) y s |x s + λ y t |x t ,

(20) which corresponds to the use of probabilistic labels. See Figure 8 in the appendix for a graphical representation of the proposed homotopy deformation. The corresponding label vector for the deformed image represented in Figure 8b is y 0.5 = [0, 0, 0.5, 0, 0, 0.5, 0, 0, 0, 0], given that λ = 0.5 and that the sampled realizations of x s and x t , represented in Figures 8a and 8c , belong to class 2 and 5, respectively.

6. EXPERIMENTAL EVALUATION:

In this section, we present some experimental evaluations of homotopy methods when applied to solve supervised regression and classification tasks. As homotopy functions we adopt the ones discussed in Section 5.2. We empirically show that homotopy methods outperform random and warm-start initialization schemes in terms of numerical performance. In particular, when the target task is complex and/or, in the transfer-learning scenario, when the data distributions are significantly different, continuation methods can achieve significant speed-up compared to random and warm-start initializations. We believe that their superior numerical performance relies on the use of homotopy functions that progressively deform the data distribution from an easy-to-solve or alreadysolved task to the target data distribution. In addition, consistently across all the benchmarks, our homotopy-based method shows faster convergence than random-initialization and faster or comparable convergence than warm-start initialization. When the source task is "similar" to the target one, there is indeed no need to gradually vary the λ parameter in Algorithm 1, but it suffices to directly set it to 1. In this extreme case, our homotopy method boils down to warm-start initialization.

6.1. REGRESSION:

For the supervised regression scenario, the problem we address is how to transfer "optimality knowledge" across two tasks that involve regressing from the input to the output of two sine wave functions with different values of phase ω. Each considered dataset has 10000 samples split across training and testing, where x and y are defined as follows

x ∼ U(0, 1) , y = sin(ωx) + ε , ε ∼ N (0, 0.01) .

The goal is to start with an "easy-to-learn" task, i.e. ω ≈ 1 rad, whose optimum is available by performing only few epochs with a first-order optimizer, e.g. SGD, Adam, and progressively transfer the "optimality knowledge" to a more complex task, i.e. ω >> 1 rad, by approximately solving the homotopy problems for increasing values of λ as described in Algorithm 1. We set ω = 1 rad for our source task distribution, and study the performance of the proposed approach with homotopy function as described in Equation ( 16) for different target distributions with ω >> 1 rad. See Figures 5a and 5b in the appendix for a visualization of the source data distribution with ω = 1 rad and the target data distribution when ω = 137 rad, respectively. The regressor is a feedforward neural network with 6 hidden layers of 100 units each and relu as activation function. In order to make the experiments more robust with respect to the choice of the step size α, we use Adam as optimizer. For the experiments in Figures 1a-1b , Figures 7a-7b in the appendix, and Figure 2a , we set α = 0.001, γ = 10, k = 200 and then performed an additional 500 epochs on the final target problem, while for the experiments in Figure 2b , we set γ = 10, k = 300 and performed an additional 600 epochs on the final target problem. In this last scenario we set α = 0.001 and then decrease it with a cosine annealing schedule to observe convergence to an optimum. As shown in Figures 1a-1b , Figures 7a-7b in the appendix, and Figures 2a and 2b , the homotopy method leads to faster convergence than the considered baselines by preserving the vicinity to an optimal solution for problems H(θ, λ) across the different λ values. In particular, we achieve a training loss up to two orders of magnitude better than the considered baselines. Figure 2 : Comparison of homotopy method, warm start and random initialization on sine wave regression tasks. The shaded areas represent the 25th and 75th percentiles. On the left, the median train loss achieved by the considered methods after 2500 epochs across 100 runs versus different omega values for the target task is plotted. For the homotopy method and warm-start initialization, ω s = 1 rad is used. On the right, the median train loss across 100 runs versus epochs for target task with ω = 137 rad is plotted. With respect to Figure 1b , in Figure 2b a cosine decay schedule is used for the learning rate, and more epochs are performed to better observe the convergence properties of the different methods.

6.2. CLASSIFICATION:

For the supervised classification scenario, we first apply the continuation method with the homotopy deformation described in Equations ( 19) and ( 20) in order to transfer optimality from the MNIST task, a notoriously "easy-to-learn" task for neural networks, to the FashionMNIST task. Since the two datasets have the same input dimensionality and the same number of classes, no additional preprocessing of the data is required. As network architecture, we use a VGG-type network, (Simonyan & Zisserman, 2015) , and Adam as optimizer with a step size of α = 0.001.

Secondly, we consider CIFAR-10 as target data distribution. Differently from the previous scenario, padding of the MNIST samples is required in order to apply Equation ( 19). The MNIST samples are also replicated across three channels. Also in this case we adopt a VGG-type network, (Simonyan & Zisserman, 2015) , and Adam as optimizer with a step size of α = 0.0001.

As shown in Figures 3a and 3b , in both benchmarks the homotopy method leads to faster convergence than random initialization. While in the second benchmark our method reaches a lower value of training loss in fewer epochs than warm-start, in the MNIST-to-FashionMNIST case the performance is comparable to using warm-start initialization. A possible interpretation is that, when the source and target task distributions are "too similar", as we hypothesize in the MNIST- to-FashionMNIST scenario, then there is no need for homotopy deformations to be applied, i.e. 0 < λ < 1, but we can directly apply λ = 1 in our scheme, which corresponds to simply using warm-start initialization.

7. CONCLUSIONS:

In this paper we propose a new methodology based on homotopy methods in order to transfer knowledge across different task distributions. In particular, our homotopy-based method allows one to exploit easy-to-solve or already-solved learning problems to solve new and complex tasks, by approximately and sequentially solving a sequence of optimization problems where the task distribution is gradually deformed from the source to the target one. We conduct a theoretical analysis of a general homotopy method in a simplified setting, and then we test our method on some popular deep learning benchmarks, where it shows superior numerical performance compared to random and warm-start initialization schemes. The proposed framework, in its limiting case, corresponds to the widely used fine-tuning heuristic, allowing for a new and more rigorous interpretation of the latter. Finally, the generality of homotopy methods also opens many novel and promising research directions in fundamental fields for deep learning, such as stochastic non-convex optimization and transfer-learning.

A PROPERTIES OF HOMOTOPIC FUNCTIONS Among the numerous properties of homotopic functions, we recall the following ones Proposition A.1. Suppose that there exists a homotopy H : Z × [0, 1] → Y from g to f , i.e. g f . Then

• g g (reflexive property)

• g f =⇒ f g (symmetric property)

• g f and f h =⇒ g h (transitive property)

Proof. See proof of Theorem 1.5 in (Suciu, 2016).

Proposition A.2. Let g, g : Z → Y and f, f : Y → W be continuous maps, and let f • g, f • g : Z → W be the respective composite maps. If g g and f f , then f • g f • g .

Proof. See proof of Proposition 1.7 in (Suciu, 2016).

B APPROXIMATION VIA GAUSSIAN FILTER:

For the supervised regression scenario, we propose the following homotopy deformation

EQUATION

A downside of this homotopy function is that the same support for x is required (the absence of the subscripts s and t on x stands to indicate that the same realization for x s and x t has to be considered). Alternatively, it is possible to approximate Equation ( 22) by using a Gaussian filter, as depicted in Figure 6b .

In particular, having sampled one realization z of the pair (x s , y s ) from the training set D s , 0 < M GF ≤ N realizations of the pair (x t , y t ) are sampled from D t . Each y t,j realization is then weighted based on the vicinity of x t,j to the sampled x s,z realization. This leads to the following approximation of the z realization of y λ

EQUATION

EQUATION

where ξ > 0 is the standard deviation of the Gaussian filter.

C ADDITIONAL FIGURES Figure 4 : Two different homotopy deformations between the probability density functions of two one-dimensional Gaussian distributions with mean and standard deviation given by µ 1 = 1, σ 1 = 1 and µ 2 = 5, σ 2 = 0.1, respectively. The homotopy represented in Figure 4a results in a mixture of Gaussian distributions, with mixture coefficient given by the homotopy parameter λ. In Figure 4b the deformation concerns instead the parameters µ and σ of the original distributions. Preserving unimodality is a desirable property when the homotopy function is used in combination with a continuation method since, as shown in Figure 4b , the location of the optimum moves together with the function deformation, allowing the optimizer to track it and gradually reach the optimum of the final target task. On the contrary, deforming the function as shown in Figure 4a does not lead to a gradual shift of the optimal solutions. Consequently, approximately and sequentially solving the problems corresponding to intermediate values of the homotopy parameter λ, i.e. 0 < λ < 1, will not allow the homotopy method to gradually approach the desired final optimal solution. (a) Homotopy transformation described in Equation ( 16).

(b) Approximation of the homotopy transformation in Equation ( 16) (also Equation ( 22)) with a Gaussian filter as described in Equations ( 23) and ( 24). Figure 8 : Graphical representation of the homotopy transformation from x s to x t as described in Equation ( 19) for two sampled images from the MNIST and FashionMNIST datasets.

D LOCAL ERROR BOUNDS FOR SGD ITERATES:

Before proving local error bounds for SGD iterates in the considered framework, given the local nature of our assumptions, we need to demonstrate two important facts, on which the proof relies.

In particular, we need to show:

• local linear contraction of Gradient Descent (GD) iterates, and that • starting in a hypersphere of radius B around a minimizer and given a "big enough" batch size, the next SGD iterate is also contained in this region for all possible realizations of the gradient estimate.

Considering problem (4) with fixed parameter λ i , in the following subsections we will refer to

θ * = θ * i , θ k = θ i,k and g k = g(θ k , λ i )

, where we drop the subscript i and the explicit dependence on λ i in order to simplify the notation. The analysis holds for all fixed parameters λ i .

D.1 LOCAL LINEAR CONTRACTION OF GD ITERATES:

Let us use GD to solve the following optimization problem

θ * ∈ arg min θ H(θ, λ i ) ,

where the objective function H fulfills Assumptions 4.2 and 4.3.

We now derive error bounds on the iterates of GD

θ k+1 = θ k -α∇ θ H(θ k , λ i ) ,

where θ k ∈ B B,θ * and 0 < α ≤ 1 L is the step size. We start by applying the definition of GD iterates and then we exploit the introduced assumptions

θ k+1 -θ * 2 = θ k -α∇ θ H(θ k , λ i ) -θ * 2 = θ k -θ * 2 -2α∇ θ H(θ k , λ i ) T (θ k -θ * ) + α 2 ∇ θ H(θ k , λ i ) 2 strong convexity ≤ (1 -αµ) θ k -θ * 2 -2α(H(θ k , λ i ) -H(θ * , λ i )) + α 2 ∇ θ H(θ k , λ i ) 2 corollary 4.2.1 ≤ (1 -αµ) θ k -θ * 2 -2α(1 -αL)(H(θ k , λ i ) -H(θ * , λ i )) .

Since H(θ k , λ i ) -H(θ * , λ i ) ≥ 0 and -2α(1 -αL) ≤ 0 when 0 < α ≤ 1 L , we can safely drop the second term and obtain the final result

θ k+1 -θ * 2 ≤ (1 -αµ) θ k -θ * 2 .

See also Theorem 2.3 in (Gower, 2018) for a derivation where Assumptions 4.2 and 4.3 are required to hold globally.

D.2 REALIZATION OF THE SGD ITERATES IN THE STRONG CONVEXITY AND L-SMOOTHNESS REGION AROUND A MINIMIZER:

We address the following optimization problem

θ * ∈ arg min θ 1 N N j=1 j (θ, λ i ) :=H(θ,λi)

, where H fulfills Assumptions 4.2-4.4.

As proved in Section D.1, under Assumptions 4.2 and 4.3, whenever θ 0 ∈ B B,θ * and 0 < α ≤ 1 L , deterministic gradient descent iterates converge linearly with contraction rate κ d := (1 -αµ).

In particular, the following inequality holds

θ D k+1 -θ * ≤ κ d • θ k -θ * ,

for any θ k such that θ k -θ * ≤ B, and superscript D denotes iterates obtained by applying the full gradient ∇H

k := ∇H(θ k , λ i ) θ D k+1 = θ k -α∇H k .

Let θ k+1 denote the iterate obtained by applying one iteration of stochastic gradient descent

θ k+1 = θ k -αg k ,

where

g k := 1 M j∈M ∇ j (θ k , λ i ) and M is a set of 0 < M ≤ N indexes randomly sampled from N = {1, . . . , N }.

Given any realization of θ k s.t. θ k -θ * ≤ B and any realization of g k , by exploiting Assumption 4.4 and the results derived in Section D.1, we have that

EQUATION

Since we have assumed that the current realization of θ k lies in the hypersphere of radius B around the optimal solution θ * , by solving for N -M N the following inequality

κ d B + 2α (N -M ) N ν ≤ B ,

we obtain that, whenever (N -M ) N ≤ (1-κ d ) 2αν B, the realization of θ k+1 will also lie in this region. These derivations show that when the realization of the current iterate θ k lies in the hypersphere of radius B around the minimizer θ * , and

(N -M ) N ≤ (1-κ d )

2αν B, then the next iterate θ k+1 will also lie in this region. Consequently, in our scenario, if we assume that the initial point θ 0 lies in the hypersphere of radius B around the minimizer θ * , then, by applying the derivations recursively, we can show that the iterates will remain in this local region around the minimizer where strong convexity and smoothness hold. where the objective function H fulfills Assumptions 4.2-4.4. We now derive error bounds for the iterates of SGD

θ k+1 = θ k -αg k ,

where g k is the unbiased estimate of ∇H k defined in the previous section and fulfills Assumption 4.5, θ k ∈ B B,θ * , 0 < α ≤ min 1 2µ , 1 L is the step size and the batch size is set to a value M such that

(N -M ) N ≤ (1-κ d ) 2αν B.

We start by applying the definition of SGD iterates

θ k+1 -θ * 2 SGD iterate = θ k -αg k -θ * 2 = θ k -θ * 2 -2αg T k (θ k -θ * ) + α 2 g k 2 .

We now take the expectation w.r.t. θ 0 , g 0 , . . . , g k-1 , g k and, considering Assumptions 4.2-4.5, we obtain the following series of inequalities

E θ0,g0,...,g k-1 ,g k θ k+1 -θ * 2 = E θ0,g0,...,g k-1 ,g k θ k -θ * 2 -2αg T k (θ k -θ * ) +α 2 g k 2 law of iterated expectations = E θ0,g0,...,g k-1 E g k θ k -θ * 2 -2αg T k (θ k -θ * ) + α 2 g k 2 | θ 0 , g 0 , . . . , g k-1 unbiased g k +bounded "variance" ≤ E θ0,g0,...,g k-1 θ k -θ * 2 -2α∇H T k (θ k -θ * ) + α 2 C 2 strong convexity ≤ (1 -2αµ) • E θ0,g0,...,g k-1 θ k -θ * 2 + α 2 C 2 .

By applying this result recursively, we derive the following bound on the error for the SGD iterates

E θ0,g0,...,g k-1 ,g k θ k+1 -θ * 2 ≤ (1 -2αµ) k+1 • E θ0 θ 0 -θ * 2 + αC 2 2µ .

See also Section 3 in (Schmidt, 2014) 

EQUATION

≤ 1 δ (1 -κ k d )B - (N -M ) N 2αν (1 -κ d )

.

Then, after applying k iterations of SGD, we obtain that

θ i+1 -θ * i+1 ≤ B -δ .

Proof.

θ i -θ * i+1 = θ i -θ * i + θ * i -θ * i+1 Triangle Ineq. ≤ θ i -θ * i + θ * i -θ * i+1

Assumption 4.7 ≤ θ i -θ * i + δ|λ i -λ i+1 | . Finally, using the fact that |λ i -λ i+1 | ≤ , it follows that, if θ i -θ * i ≤ B -δ with 0 ≤ ≤ B δ , then θ i -θ * i+1 ≤ B. We now derive the conditions on such that θ i+1 -θ * i+1 ≤ B -δ . By applying recursively the results derived in Section D.2 (25), we obtain that

θ i+1 -θ * i+1 ≤ κ k d θ i -θ * i+1 + 2α (N -M ) N ν k-1 i=0 κ i d .

By using the limit of the geometric series, we have that

θ i+1 -θ * i+1 ≤ κ k d θ i -θ * i+1 + (N -M ) N 2αν (1 -κ d )

.

Finally, by considering that θ i -θ * i+1 ≤ B and by solving in the following inequality

κ k d B + (N -M ) N 2αν (1 -κ d ) ≤ B -δ ,

we obtain the following upper bound on

≤ 1 δ (1 -κ k d )B - (N -M ) N 2αν (1 -κ d ) ,

from which also the extra condition on the batch size (28)

(N -M ) N ≤ (1 -κ k d )(1 -κ d )B 2αν .

Proof.

E θ i+1 -θ * i+1 2 Ineq. 10 ≤ (1 -2αµ) k E θ i -θ * i+1 2 + αC 2 2µ = (1 -2αµ) k E θ i -θ * i + θ * i -θ * i+1 2 + αC 2 2µ Triangle Ineq. ≤ (1 -2αµ) k E θ i -θ * i + θ * i -θ * i+1 2 + αC 2 2µ = (1 -2αµ) k E θ i -θ * i 2 + θ * i -θ * i+1 2 +2 θ i -θ * i θ * i -θ * i+1 + αC 2 2µ Assumption 4.7 ≤ (1 -2αµ) k E θ i -θ * i 2 + δ 2 |λ i -λ i+1 | 2 +2δ θ i -θ * i |λ i -λ i+1 |)] + αC 2 2µ ≤ (1 -2αµ) k δ 2 ˜ 2 + 2δr θ ˜ + r 2 θ + αC 2 2µ .

We now solve in ˜ the following second degree inequality

EQUATION

The inequality (29) admits solutions if and only if r 2 θ ≥ αC 2 2µ . In particular, inequality (29) holds ∀˜ ∈ [0, ¯ ], where ¯ = -r θ δ + 1 δ r 2 θ -αC 2 /2µ

(1-2αµ) k .

F EXPERIMENTAL EVALUATION: TEST PERFORMANCES 

:



Figures and Tables:

(a) ωs = 1 rad, ωt = 74 rad.(b) ωs = 1 rad, ωt = 137 rad.

Figure 1: Median train loss across 100 runs versus epochs for sine wave regression tasks with different frequency values. The shaded areas represent the 25th and 75th percentiles. See Section F.1 in the appendix for an evaluation of the test performance.

Figure3: Median train loss across 10 runs versus epochs for different target task distributions. In both cases, the source task is the classification of the MNIST dataset. See Section F.2 in the appendix for an evaluation of the test performance.

(a) ω = 1 rad. (b) ω = 137 rad.

Figure 5: Graphical representations of the source (left) and target with ω = 137 rad (right) data distributions used for the sine-wave regression evaluation.

Figure 6: Graphical representation of the proposed homotopy transformation for the supervised regression scenario when applied to progressively deform a sine wave function with frequency of 1 radian into a sine wave function with frequency of 137 radians for different values of homotopy parameter.

Figure 7: Median train loss across 100 runs versus epochs for sine wave regression tasks with different omega values.

PROOF OF PROPOSITION 4.9 Let us use SGD to solve the following optimization problem θ * ∈ arg min θ H(θ, λ i ) ,

Figure 9: Graphical representation of the results derived in Proposition 4.11. The continuous and dashed lines are used to represent the circles of radius B and B -δ around the optimal solutions, respectively.

ωs = 1 rad, ωt = 74 rad. (b) ωs = 1 rad, ωt = 95 rad. (c) ωs = 1 rad, ωt = 116 rad. (d) ωs = 1 rad, ωt = 137 rad.

Figure 10: Median test loss across 100 runs versus epochs for target tasks with different ω values.The shaded areas represent the 25th and 75th percentiles. For warm-start initialization and homotopy method, ω s = 1 rad is used for the source task.

for a derivation where Assumptions 4.2 and 4.3 are required to hold globally.

Table content: <table/>



Table content: <table><tr><td>F.2 CLASSIFICATION</td><td/><td/></tr><tr><td>Method</td><td colspan="2">Final Mean Test Accuracy Best Mean Test Accuracy</td></tr><tr><td>homotopy γ = 5, k = 2</td><td>0.89 ± 0.003</td><td>0.91 ± 0.002</td></tr><tr><td>homotopy γ = 5, k = 4</td><td>0.89 ± 0.002</td><td>0.91 ± 0.003</td></tr><tr><td>homotopy γ = 10, k = 1</td><td>0.89 ± 0.004</td><td>0.91 ± 0.001</td></tr><tr><td>homotopy γ = 10, k = 4</td><td>0.90 ± 0.002</td><td>0.91 ± 0.003</td></tr><tr><td>warm start</td><td>0.89 ± 0.003</td><td>0.90 ± 0.002</td></tr><tr><td>random init</td><td>0.89 ± 0.004</td><td>0.90 ± 0.003</td></tr><tr><td/><td>Table 1: MNIST-FashionMNIST</td><td/></tr><tr><td>Method</td><td colspan="2">Final Mean Test Accuracy Best Mean Test Accuracy</td></tr><tr><td>homotopy γ = 5, k = 2</td><td>0.55 ± 0.004</td><td>0.59 ± 0.003</td></tr><tr><td>homotopy γ = 10, k = 1</td><td>0.55 ± 0.005</td><td>0.60 ± 0.002</td></tr><tr><td>homotopy γ = 10, k = 2</td><td>0.56 ± 0.003</td><td>0.60 ± 0.003</td></tr><tr><td>homotopy γ = 10, k = 4</td><td>0.56 ± 0.005</td><td>0.61 ± 0.004</td></tr><tr><td>warm start</td><td>0.54 ± 0.006</td><td>0.59 ± 0.005</td></tr><tr><td>random init</td><td>0.64 ± 0.02</td><td>0.64 ± 0.02</td></tr></table>

Peer review (iclr2020):
Based on homotopy,, the paper describes a more rigorous approach to transfer learning than the so called ‘fine-tuning’ heuristic. Progress in the direction of more principled approaches for transfer learning would be tremendously impactful, since one of the core promises of deep learning is the learning of features, which can be used in different downstream tasks. Essentially, (if this reviewer understood this correctly) the idea behind this paper works by interpolation between the original task of interest and a potentially easier to optimize surrogate task. Overall, this reviewer found the concept simple and elegant, and well motivated, and also well introduced. However, since this reviewer does not have a formal background in mathematics, they cannot assess the soundness of the proofs. The paper tests the hypothesis by a simple function approximation regression task, and a classification task to learn to transfer from MNIST to fashion MNIST and MNIST to CIFAR, with promising results. One might argue that a more thorough evaluation would have been desirable, since the claims made by the paper are quite general, and it would have been in the authors’ best interest to present more thorough evidence that their concept works on wider scale of problems, ideally on an NLP task, given the current hype on pre-training with Transformer-based models. Previous work & citations: I would recommend to cite Schmidhuber 1987 (Evolutionary principles in self-referential learning) and Hochreiter et al 2001 (Learning to Learn with gradient descent) in the context of Meta learning. It would be nice to cite Klambauer et al (Self normalizing Networks) in the context of speeding up deep neural network training. The citations of the VGG paper is currently referenced by first names of the authors, not their last names, I am not sure if this was intended.