TRANSFERRING OPTIMALITY ACROSS DATA DISTRI-
Year: 2020

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.

heading: 1 INTRODUCTION
text: 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).

heading: 2 RELATED WORK
text: 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), Kä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.

heading: 3 BACKGROUND
text: In this work, we will focus on solving problems of the form
θ∗ ∈ arg min θ∈Rd
1
N N∑ j=1
`j(θ)︸ ︷︷ ︸ :=J(θ) , (1)
where J : Rd → 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.

heading: 3.1 HOMOTOPIC FUNCTIONS AND CONTINUATION METHODS FOR OPTIMIZATION
text: 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 H : Z × [0, 1]→ Y such that
H(z, 0) = g(z) , H(z, 1) = f(z) , ∀z ∈ Z . (2)

heading: If such function H exists, g is said to be homotopic of f , and this relation is denoted by g ' f .
text: It is straightforward to show that, A ⊆ Rn 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) = J0(θ) 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
θ∗i ∈ arg min θ∈Rd
1
N N∑ j=1
`j(θ, λi)︸ ︷︷ ︸ :=H(θ,λi) , (3)
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.

heading: 4 THEORETICAL ANALYSIS
text: 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 θ∈Rd H(θ, λi) , ∀i = 1, . . . , γ , (4)
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.

heading: 4.1 ASSUMPTIONS
text: 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 BB,θ∗i be the following set of vectors
BB,θ∗i := {θ s.t. ‖θ − θ ∗ i ‖ ≤ B} , ∀i = 0, . . . , γ .
Assumption 4.2 (local L-smoothness). Assume that there exists a constant L > 0 such that
‖∇θH(θ̃, λi)−∇θH(θ̂, λi)‖ ≤ L‖θ̃ − θ̂‖ , ∀θ̃, θ̂ ∈ BB,θ∗i , ∀i = 0, . . . , γ . (5)
Corollary 4.2.1. If H is locally L-smooth in θ, then the following inequality holds
H(θ∗i , λi)−H(θ̂, λi) ≤ − 1
2L ‖∇θH(θ̂, λi)‖2 , ∀θ̂ ∈ BB,θ∗i , ∀i = 0, . . . , γ. (6)
Proof. See Lemma 1.1 in (Gower, 2018) for a proof.
Assumption 4.3 (local µ-strong convexity). Assume that there exists µ > 0 such that
H(θ̃, λi) ≥ H(θ̂, λi) +∇θH(θ̂, λi)T (θ̃− θ̂) + µ
2 ‖θ̃− θ̂‖2 , ∀θ̃, θ̂ ∈ BB,θ∗i , ∀i = 0, . . . , γ . (7)
Assumption 4.4 (bounded `-derivative). Assume that there exists ν > 0 such that
‖∇θ`j(θ̂, λi)‖ ≤ ν , ∀θ̂ ∈ BB,θ∗i , ∀i = 0, . . . , γ, ∀j = 1, . . . , N. (8)
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
E [ ‖g(θ̂, λi)‖2 ] ≤ C2 , ∀θ̂ ∈ BB,θ∗i , ∀i = 0, . . . , γ. (9)
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)).

heading: 4.2 FUNDAMENTAL THEORETICAL PRELIMINARIES
text: 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 ∈ BB,θ∗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
E [ ‖θi+1 − θ∗i+1‖2 ] ≤ (1− 2αµ)k · E [ ‖θi − θ∗i+1‖2 ] + αC2
2µ . (10)
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.

heading: 4.3 MAIN THEORETICAL CONTRIBUTIONS
text: 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 ] ≤ r2θ , 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 ] ≤ r2θ . 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 ∈ BB,θ∗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 ∈ BB,θ∗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− κkd)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 κd = √
(1− αµ). For θ0 ∈ BB−δ ,θ∗0 and rθ ∈ R such that
r2θ ≥ αC2
2µ , (11)
then, if E [ ‖θi − θ∗i ‖2 ] ≤ r2θ and |λi − λi+1| ≤ ̃, where ̃ := min {̄, } with
̄ = −rθ δ + 1 δ √ r2θ − αC2/2µ (1− 2αµ)k , (12)
the following inequality holds E [ ‖θi+1 − θ∗i+1‖2 ] ≤ r2θ . (13)
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. r2θ ≥ α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.

heading: 5 TRANSFERRING OPTIMALITY VIA HOMOTOPY METHODS
text: 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.

heading: 5.1 PROBLEM FORMULATION
text: 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 θ ∈ Rd. In the supervised learning scenario, independently from the type of task t, we typically dispose of a training set Dt consisting of N pairs of examples (xj , yj). 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 ŷ = f(x; θ) and the target outputs y. In particular, the learning process consists in minimizing the following empirical objective function
J(θ) := 1
N ∑ (xj ,yj)∈Dt `(yj , f(xj ; θ)) , (14)
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 Qx,y , with density qx,y = q(y|x)q(x), and the learned data distribution Px,y(θ), with density px,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
J(θ) = 1
N ∑ (xj ,yj)∈Dt q(y|x) log q(y|x) p(y|x; θ) . (15)

heading: 5.2 HOMOTOPY FUNCTIONS ACROSS DATA DISTRIBUTIONS
text: 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 Qx,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 Ds of pairs (xs, ys) ∼ Qxs,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 Dt of pairs (xt, yt) ∼ Qxt,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.

heading: 5.2.1 SUPERVISED REGRESSION
text: 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
yλ|x = (1− λ) ys|x+ λ yt|x , (16) p(yλ|x) = N (yλ ; f(x; θ), σ2 I) . (17)
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
H(θ, λ) := E(x,yλ) ‖(1− λ) (ys − f(x; θ)) + λ (yt − f(x; θ)) ‖ 2 . (18)
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 xs and xt 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).

heading: 5.2.2 SUPERVISED CLASSIFICATION
text: In the case of supervised classification, by modeling the density of the conditional learned distribution as p(y|x; θ) = Multinoulli(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
xλ = (1− λ)xs + λxt , (19) yλ|xλ = (1− λ) ys|xs + λ yt|xt , (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 y0.5 = [0, 0, 0.5, 0, 0, 0.5, 0, 0, 0, 0], given that λ = 0.5 and that the sampled realizations of xs and xt, represented in Figures 8a and 8c, belong to class 2 and 5, respectively.

heading: 6 EXPERIMENTAL EVALUATION
text: 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.

heading: 6.1 REGRESSION
text: 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) . (21)
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.

heading: 6.2 CLASSIFICATION
text: 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.

heading: 7 CONCLUSIONS
text: 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.

heading: ACKNOWLEDGMENTS
text: This work has partly been supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant no. 716721 as well as by the German Federal Ministry for Economic Affairs and Energy (BMWi) via DyConPV (0324166B), and by DFG via Research Unit FOR 2401. In addition, Q. Tran-Dinh has partly been supported by the National Science Foundation (NSF), grant. no. 1619884. The authors thank Stefan Falkner for his helpful suggestions and comments.

heading: A PROPERTIES OF HOMOTOPIC FUNCTIONS
text: 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).

heading: B APPROXIMATION VIA GAUSSIAN FILTER
text: For the supervised regression scenario, we propose the following homotopy deformation
yλ|x = λ ys|x+ (1− λ) yt|x . (22)
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 xs and xt 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 (xs, ys) from the training set Ds, 0 < MGF ≤ N realizations of the pair (xt, yt) are sampled from Dt. Each yt,j realization is then weighted based on the vicinity of xt,j to the sampled xs,z realization. This leads to the following approximation of the z realization of yλ
yλ,z = (1− λ) ys,z + λ
MGF MGF∑ j=1 wj yt,j , (23)
wj = 1√ 2πξ2 exp
( −||xs,z − xt,j || 2
2ξ2
) , (24)
where ξ > 0 is the standard deviation of the Gaussian filter.

heading: C ADDITIONAL FIGURES

heading: D LOCAL ERROR BOUNDS FOR SGD ITERATES
text: 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 gk = 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 ∈ BB,θ∗ and 0 < α ≤ 1L 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 < α ≤ 1L , 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 ∈ BB,θ∗ and 0 < α ≤ 1L , deterministic gradient descent iterates converge linearly with contraction rate κd := √ (1− αµ).
In particular, the following inequality holds
‖θDk+1 − θ∗‖ ≤ κd · ‖θk − θ∗‖ ,
for any θk such that ‖θk− θ∗‖ ≤ B, and superscript D denotes iterates obtained by applying the full gradient∇Hk := ∇H(θk, λi)
θDk+1 = θk − α∇Hk .
Let θk+1 denote the iterate obtained by applying one iteration of stochastic gradient descent
θk+1 = θk − αgk , where gk := 1M ∑ 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 gk, by exploiting Assumption 4.4 and the results derived in Section D.1, we have that
‖θk+1 − θ∗‖ = ‖θk − αgk − θ∗‖ = ‖θk − α∇Hk + α∇Hk − αgk − θ∗‖ ≤ ‖θk − α∇Hk − θ∗‖+ α‖∇Hk − gk‖
= ‖θk − α∇Hk − θ∗‖+ α ∥∥∥ 1 N ∑ j∈N\M ∇`j + 1 N ∑ j∈M ∇`j − 1 M ∑ j∈M ∇`j ∥∥∥
= ‖θk − α∇Hk − θ∗‖+ α ∥∥∥ 1 N ∑ j∈N\M ∇`j + M −N NM ∑ j∈M ∇`j ∥∥∥
≤ ‖θk − α∇Hk − θ∗‖+ α  1 N ∑ j∈N\M ‖∇`j‖+ N −M NM ∑ j∈M ‖∇`j‖  ≤ ‖θDk+1 − θ∗‖+ 2α
(N −M) N ν
≤ κd‖θk − θ∗‖+ 2α (N −M)
N ν .
(25)
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−MN the following inequality
κdB + 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.
D.3 PROOF OF PROPOSITION 4.9
Let us use SGD to solve the following optimization problem
θ∗ ∈ arg min θ H(θ, λi) ,
where the objective function H fulfills Assumptions 4.2– 4.4. We now derive error bounds for the iterates of SGD θk+1 = θk − αgk , where gk is the unbiased estimate of ∇Hk defined in the previous section and fulfills Assumption 4.5, θk ∈ BB,θ∗ , 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 − αgk − θ∗‖2
= ‖θk − θ∗‖2 − 2αgTk (θk − θ∗) + α2‖gk‖2 . We now take the expectation w.r.t. θ0, g0, . . . , gk−1, gk and, considering Assumptions 4.2- 4.5, we obtain the following series of inequalities
Eθ0,g0,...,gk−1,gk [ ‖θk+1 − θ∗‖2 ] = Eθ0,g0,...,gk−1,gk [ ‖θk − θ∗‖2 − 2αgTk (θk − θ∗)
+α2‖gk‖2 ]
law of iterated expectations = Eθ0,g0,...,gk−1 [ Egk [ ‖θk − θ∗‖2
−2αgTk (θk − θ∗) + α2‖gk‖2 | θ0, g0, . . . , gk−1 ]]
unbiased gk+bounded “variance” ≤ Eθ0,g0,...,gk−1 [ ‖θk − θ∗‖2
−2α∇HTk (θk − θ∗) ] + α2C2
strong convexity ≤ (1− 2αµ) · Eθ0,g0,...,gk−1 [ ‖θk − θ∗‖2 ] + α2C2 .
By applying this result recursively, we derive the following bound on the error for the SGD iterates
Eθ0,g0,...,gk−1,gk [ ‖θk+1 − θ∗‖2 ] ≤ (1− 2αµ)k+1 · Eθ0 [ ‖θ0 − θ∗‖2 ] + αC2
2µ .
See also Section 3 in (Schmidt, 2014) for a derivation where Assumptions 4.2 and 4.3 are required to hold globally.

heading: E MAIN THEORETICAL CONTRIBUTIONS
text: E.1 PROOF OF PROPOSITION 4.11
Proposition E.1. Let θi ∈ BB,θ∗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− κkd)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‖ ≤ κkd‖θi − θ∗i+1‖+ 2α (N −M)
N ν k−1∑ i=0 κid .
By using the limit of the geometric series, we have that
‖θi+1 − θ∗i+1‖ ≤ κkd‖θi − θ∗i+1‖+ (N −M)
N
2αν
(1− κd) .
Finally, by considering that ‖θi − θ∗i+1‖ ≤ B and by solving in the following inequality
κkdB + (N −M)
N
2αν
(1− κd) ≤ B − δ ,
we obtain the following upper bound on
≤ 1 δ
( (1− κkd)B −
(N −M) N 2αν (1− κd)
) ,
from which also the extra condition on the batch size (N −M)
N ≤ (1− κ k d)(1− κd)B 2αν .
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.
E.2 PROOF OF THEOREM 4.12
Theorem E.2. 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 κd = √
(1− αµ). For θ0 ∈ BB−δ ,θ∗0 and rθ ∈ R such that
r2θ ≥ αC2
2µ , (26) then, if E [ ‖θi − θ∗i ‖2 ] ≤ r2θ and |λi − λi+1| ≤ ̃, where ̃ := min {̄, } with
̄ = −rθ δ + 1 δ √ r2θ − αC2/2µ (1− 2αµ)k , (27)
the following inequality holds E [ ‖θi+1 − θ∗i+1‖2 ] ≤ r2θ . (28)
Proof.
E [ ‖θi+1 − θ∗i+1‖2 ] Ineq. 10 ≤ (1− 2αµ)kE [ ‖θi − θ∗i+1‖2 ] + αC2
2µ = (1− 2αµ)kE [ ‖θi − θ∗i + θ∗i − θ∗i+1‖2 ] + αC2
2µ Triangle Ineq. ≤ (1− 2αµ)kE [( ‖θi − θ∗i ‖+ ‖θ∗i − θ∗i+1‖ )2] + αC2
2µ = (1− 2αµ)kE [( ‖θi − θ∗i ‖2 + ‖θ∗i − θ∗i+1‖2
+2‖θi − θ∗i ‖‖θ∗i − θ∗i+1‖ )] + αC2
2µ Assumption 4.7 ≤ (1− 2αµ)kE [( ‖θi − θ∗i ‖2 + δ2|λi − λi+1|2
+2δ‖θi − θ∗i ‖|λi − λi+1|)] + αC2
2µ ≤ (1− 2αµ)k ( δ2̃2 + 2δrθ ̃+ r 2 θ ) + αC2
2µ .
We now solve in ̃ the following second degree inequality
(1− 2αµ)k ( δ2̃2 + 2δrθ ̃+ r 2 θ ) + αC2
2µ ≤ r2θ . (29)
The inequality (29) admits solutions if and only if r2θ ≥ αC 2 2µ . In particular, inequality (29) holds
∀̃ ∈ [0, ̄], where ̄ = − rθδ + 1 δ √ r2θ−αC2/2µ (1−2αµ)k .

heading: F EXPERIMENTAL EVALUATION: TEST PERFORMANCES
text: F.1 REGRESSION
F.2 CLASSIFICATION

Review summaries:
This paper proposes algorithm for transferring knowledge from easy -to-solved to complex tasks or from already solve to new tasks. It relies on homotopy functions and sequentially solves a sequence of optimization problems where the task distribution is gradually deformed from a source task to the target task. Theoretical guarantees are provided and proven in a strongly convex setting. The main results from the theory show that the distance between the final solution and its optimal are less or equal to  relative to the distance of the initial source solution to its optimum. So a near optimal solution for the source task will lead to near optimal solution for the target task. Regression and Classification experimentations show competitive results compared to  random and warm-start initialization schemes.