Year: 2019

Abstract:
We present a Homotopy Training Algorithm (HTA) to solve optimization problems arising from neural networks. The HTA starts with several decoupled systems with low dimensional structure and tracks the solution to the high dimensional coupled system. The decoupled systems are easy to solve due to the low dimensionality but can be connected to the original system via a continuous homotopy path guided by the HTA. We have proved the convergence of HTA for the non-convex case and existence of the homotopy solution path for the convex case. The HTA has provided a better accuracy on several examples including VGG models on CIFAR-10. Moreover, the HTA would be combined with the dropout technique to provide an alternative way to train the neural networks.

heading: 1 INTRODUCTION
text: Deep learning methods, as the rising star among all machine learning methods in recent years, have already had great success in many applications. Many advancements Cang et al. (2018); Cang & Wei (2018); Yin et al. (2018; 2016) on deep learning have been made in the last few years. However, as the size of new state-of-the-art models keep growing larger and larger, they rely more and more heavily on efficient algorithms for training and making inferences from such models. This clearly puts strong limitations on the application scenarios of nerual netowrks on roboticsKonda et al. (2012), auto-pilot automobiles Hammerla et al. (2016), and aerial systems Maire et al. (2014). At present, there are two big challenges in fundamentally understanding deep neural networks:
• How to efficiently solve the highly nonlinear and non-convex optimization problems that arise during the training of a deep learning model?
• How to design a deep neural network structure for specific problems?
In order to solve these challenges, in this paper, we will present a new training algorithm based on homotopy continuation method Bates et al. (2011; 2008); Morgan & Sommese (1987) which has been successfully used to study nonlinear problems such as nonlinear differential equations and optimizations Hao et al. (2012); Hao & Harlim. In order to tackle the nonlinear optimization problem in neural networks, the homotopy training algorithm is designed and shows efficiency and feasibility for both convolutional and fully connected neural networks with complex structures. The idea of the HTA is to start with a simplified optimization problem and construct a continuous path to the original problem. Specifically, let’s consider the following general optimization problem
min θ f(θ) where θ ∈ Rn. (1)
We then divide θ in two groups, namely, θ = θ1∪ θ2, where θ1 ∈ Rn1 and θ2 ∈ Rn2 (n1 +n2 = n). (Although that it can be easily extended to any m decompositions, we use m = 2 as an illustration of the HTA idea.) Then we first solve the following decoupled simplified optimization problems
min θ1 f(θ1, θ̃2) and min θ2 f(θ̃1, θ2) (2)
where θ̃ = θ̃1 ∪ θ̃2 is a random chosen point in Rn for instance we choose θ̃ = 0 in the dropout Srivastava et al. (2014). Second we define the following homotopy function
H(Θ, L) = f(θ) + L‖θ1 − θ̃1‖2 + f(η) + L‖η2 − θ̃2‖2 + 1
L ‖θ − η‖2, (3)
where Θ = (θ, η) and L ∈ (0,∞) is the homotopy parameter. When L → ∞, we have two decoupled optimization problems which we have already solved in the first step; when L → 0, we recover the original optimization by forcing θ = η. By gradually increasing parameter L, we are able to construct a continuous path Θ(L)between two decoupled low-dimensional optimization problems and the original optimization problem.
Remark The HTA can be reviewed in point view of nonlinear equations by considering the necessary condition of the optimization problem (1):
H(Θ, L) = (∇θf(θ) η − θ ) 1 L + L  ∇θ2f(θ) θ1 − θ̃1 ∇η1f(η) η2 − θ̃2  . When L is large, θ2 and η1 can be solved by ∇θ2f(θ) = 0 and ∇η1f(η) = 0 respectively. When L is small, we have the solution of ∇θf(θ) = 0 by the diagonal homotopy θ − η = 0. This technique has been widely used in solving systems of nonlinear equations especially for polynomial systems Bates et al. (2008); Wampler & Sommese (2005). Since the θ̃ is randomly chosen, the universal approximation theory Park & Sandberg (1991) can guarantee that the random point θ̃ in Rn is always connected to the solution of ∇θf(θ) = 0. In another word, the solution path Θ(L) always exists or the probability that Θ(L) is not connected is zero Wampler & Sommese (2005).

heading: 2 APPLICATIONS OF HTA ON NEURAL NETWORK

heading: 2.1 CONVOLUTIONAL NEURAL NETWORKS
text: In this subsection, we will specify the variable θ and the objective function f(θ) of (3) in the convolutional neural network (CNN) setup. Definition 2.1. Suppose we have two tensors, Al×m×n = (ai,j,k)1≤i≤l,1≤j≤m,1≤k≤n, Bc×l×m′×n′ = (bh,i,j,k)1≤h≤c,1≤i≤l,1≤j≤m′,1≤k≤n′ , where c,m, n,m′, n′, l ∈ N and m′ ≤ m,n′ ≤ n. The convolution of two tensors is defined by the following:
A ∗Bh,j,k = l∑
i′=1 m′∑ j′=1 n′∑ k′=1 ai′,j′+j−1,k′+k−1 · bh,i′,j′,k′ , (4)
where B is called the convolutional kernel.
Then with the input x and the output y, the CNN is written as
y(x; θ) = σ(σ(· · ·σ(σ(x ∗W1 + b1) ∗W2 + b2) · · ·+ bn−1) ∗Wn + bn) (5) where Wi ∈ Rci,li,mi,ni and σ is the activation function. Then we have θ = {Wi, bi}ni=1 which can be easily divided into two subgroups θj = {Wij , bi}ni=1 where Wij ∈ R ci 2 ,li,mi,ni . The objective function f(θ) in (3) is defined by the cross entropy loss function Goodfellow et al. (2016) and defined as L(x, y) = −x(y) + log( ∑ j e x[j]). Due to the random algorithms in machine learning Bottou (2010); Kingma & Ba (2014), the homotopy function in (3) can be rewritten as H(Θ; ξ, L) where ξ is a random variable due to random algorithms.

heading: 2.2 FULLY CONNECTED NEURAL NETWORKS
text: Similar to the CNN, a fully connected neural network with n hidden layers is written as
y(x, θ) = σ(Wnσ(Wn−1σ(W2σ(W1x+ b1) + b2) · · ·+ bn−1) + bn), (6) where Wi ∈ Rn×di , βi ∈ Rdi , n is the dimension of the input x and di is the width of the i-th layer. Then we have θ = {Wi, bi}ni=1 which can be easily divided into two subgroups θ1 = {Wi, bi} n/2 i=1 and θ2 = {Wi, bi}ni=n/2+1. Or we can use the dropout technique to randomly split θ as θ1 and θ2. Similarly, the homotopy function in (3) can be also rewritten as H(Θ; ξ, L) due to the random algorithms Bottou (2010); Kingma & Ba (2014).

heading: 3 THEORETICAL ANALYSIS
text: In this section, we will provide a theoretical analysis of the HTA by using stochastic gradient scheme to solve (3) for each L. We redefine our stochastic gradient scheme for the homotopy process as
θk+1 = θk − γkG(θk; ξk, Lk), (7)
where G(θ; ξ, L) := ∇θF (θ; ξ, L), γk is the learning rate, L0 = M , Lk ↗ 0 and M is a given large number. Then we consider the expectation of the homotopy function as Eξ(H(θ; ξ, L)), for simplicity, and denote Eξ[H(θ; ξ, L)] as g(θ;L).

heading: 3.1 CONVERGENCE OF THE HTA FOR ANY GIVEN L IN THE NON-CONVEX CASE
text: Theorem 3.1. (Nonconvex Convergence) If {θk} is contained in a bounded open set, supposing that (7) is run with a step size sequence satisfying
∞∑ k=1 γk =∞ and ∞∑ k=1 γ2k <∞, (8)
then we have
E[ 1
Ak K∑ k=1 γk‖∇g(θk;L)‖22]→ 0 as K →∞ with Ak := K∑ k=1 γk. (9)
Proof. First, we prove the Lipschitz-continuous objective gradients condition Bottou et al. (2018), which means that g(θ;L) is C1 and ∇g(θ;L) is Lipschitz continuous:
• g(θ;L) is C1. Since y(xξ; θ) ∈ C1 and L(·, y) is C1, we have that H(θ; ξ, L) ∈ C1 or that ∇θH(θ; ξ, L) is continuous. Considering
∇θg(θ;L) = ∇θEξ(H(θ; ξ, L)) = Eξ(∇θH(θ; ξ, L)), (10)
we have that∇θg(θ, L) is continuous or that g(θ;L) ∈ C1.
• ∇g(θ;L) is Lipschitz continuous. Since ∇θg(θ;L) = Eξ(∇θH(θ; ξ, L)), we only need to prove Eξ(∇θL(y(xξ; θ), yξ)) is Lipschitz continuous. In fact
∇θL(y(xξ; θ), yξ) = ∇xL(y(xξ; θ), yξ) · ∇θy(xξ; θ) (11)
we will show that ∇xL(y(xξ; θ), yξ) and ∇θy(xξ; θ) are bounded and Lipschitz continuous. Since that both σ(x) = 11+e−x and σ
′(x) = σ(x)(1−σ(x)) are bounded and Lipschitz continuous, and that {θk} is in an open set (assumption of Theorem 3.1), ∇θy(xξ; θ) is bounded and Lipschitz continuous. (xξ is bounded because the size of our dataset is finite.) By differentiating L(x, y), we have
∇xL(x, y) = (−δ1y + ex1∑ j e xj , · · · ,−δny + exn∑ j e xj ), (12)
where δ is the Kronecker delta. since
∂
∂xk exi∑ j e xj =  exi ∑ j e xj−(exi )2 ( ∑ j e xj )2 k = i −exiexk ( ∑
j e xj )2
k 6= i, (13)
which implies that ∣∣ ∂ ∂xk exi∑ j e xj ∣∣ ≤ 2, we have ∇xL(·, y) is Lipschitz continuous and bounded. Therefore ∇xL(y(xξ; θ), yξ) is Lipschitz continuous and bounded. Therefore we have that ∇g(θ;L) is Lipschitz continuous.
Second, we prove the first and second moment limits condition Bottou et al. (2018):
a. According to our theorem assumption, {θk} is contained in an open set which is bounded. Since that g is continuous, g is bounded;
b. Since that G(θk; ξk, L) = ∇θF (θk; ξk, L) is continuous, we have
Eξk [G(θk; ξk, L)] = ∇θEξk [H(θk; ξk, L)] = ∇θg(θk;L). (14)
Therefore
∇g(θk;L)TEξk [G(θk; ξk, L)] = ∇g(θk;L)T ·∇g(θk;L) = ‖∇g(θk;L)‖22 ≥ u‖∇g(θk;L)‖22 (15) for 0 < u ≤ 1. On the other hand, we have
‖Eξk [G(θk; ξk, L)]‖2 = ‖∇g(θk;L)‖2 ≤ uG‖∇g(θk;L)‖22 (16)
for uG ≥ 1.
c. Since that∇H(θk; ξk, L) is bounded for a given L, we have Eξk [‖∇H(θk; ξk, L)‖22] is also bounded. Thus
Vξk [G(θk; ξk, L)] := Eξk [‖G(θk; ξk, L)‖22]−‖Eξk [G(θk; ξk, L)]‖22 ≤ Eξk [‖G(θk; ξk, L)‖22], (17) which implies that Vξk [G(θk; ξk, L)] is bounded.
We have checked assumptions 4.1 and 4.3 in Bottou et al. (2018). By theorem 4.10 in Bottou et al. (2018), with the diminishing stepsize, namely,
∞∑ k=1 γk =∞ and ∞∑ k=1 γ2k <∞, (18)
the following convergence is obtained
E[ 1
Ak K∑ k=1 γk‖∇g(θk;L)‖22]→ 0 as K →∞. (19)
3.2 EXISTENCE OF SOLUTION PATH θ(L) FOR THE CONVEX CASE
Second we explore the existence of solution path θ(L) when t varies fromM to 0 for the convex case theoretically. This is quite difficult for the non-convex case although we may find it numerically. However, theoretical analysis of the convex case would also give us some guidance on how the HTA obtains an optimal solution for L = 0.
Theorem 3.2. Assume that g(·, ·) is convex and differentiable and that ‖G(θ; ξ, L)‖ ≤ M0. Then for stochastic gradient scheme (7), we have
lim n→∞
E[g(θ̄n, L̄n)] = g(θ1∗, 0), (20)
where θ̄n = ∑n
k=0 γkθk∑n k=0 γk
and L̄n = ∑n
k=0 γkLk∑n k=0 γk .
Proof.
E[‖θk+1 − θ Lk+1 ∗ ‖2] = E[‖θk − γkG(θk; ξk, Lk)− θLk∗ ‖2]
−2E[〈θk − γkG(θk; ξk, Lk)− θLk∗ , θ Lk+1 ∗ − θLk∗ 〉] +E[‖θLk+1∗ − θLk∗ ‖2]. (21)
By defining
Ak = −2E[〈θk − γkG(θk; ξk, Lk)− θLk∗ , θ Lk+1 ∗ − θLk∗ 〉] + E[‖θ Lk+1 ∗ − θLk∗ ‖2], (22)
and assuming that ∑m k=0Ak ≤ A <∞, we obtain
E[‖θk+1 − θ Lk+1 ∗ ‖2] = E[‖θk − γkG(θk; ξk, Lk)− θLk∗ ‖2] +Ak
= E[‖θk − θLk∗ ‖2]− 2γkE[〈G(θk; ξk, Lk), θk − θLk∗ 〉] + γ2kE[‖G(θk; ξk, Lk)‖2] +Ak ≤ E[‖θk − θLk∗ ‖2]− 2γkE[〈G(θk; ξk, Lk), θk − θLk∗ 〉] + γ2kM2 +Ak
Since that
E[〈G(θk; ξk, Lk), θk − θLk∗ 〉] = Eξ0,··· ,ξk−1 [Eξk [〈G(θk; ξk, Lk), θk − θLk∗ 〉|ξ0, · · · , ξk−1]] = Eξ0,··· ,ξk−1 [〈∇g(θk;Lk), θk − θLk∗ 〉|ξ0, · · · , ξk−1] = E[〈∇g(θk;Lk), θk − θLk∗ 〉], (23)
we have
E[‖θk+1 − θ Lk+1 ∗ ‖2] ≤ E[‖θk − θLk∗ ‖2]− 2γkE[〈∇g(θk;Lk), θk − θLk∗ 〉] + γ2kM2 +Ak. (24)
Due to the convexity of g(·, Lk), namely, 〈∇g(θk, Lk), θk − θLk∗ 〉 ≥ g(θk;Lk)− g(θLk∗ ;Lk), (25)
we conclude that
E[‖θk+1 − θ Lk+1 ∗ ‖2] ≤ E[‖θk − θLk∗ ‖2]− 2γkE[g(θk;Lk)− g(θLk∗ , Lk)] + γ2kM2 +Ak, (26)
or
2γkE[g(θk;Lk)− g(θLk∗ ;Lk)] ≤ −E[‖θk+1 − θ Lk+1 ∗ ‖2 − ‖θk − θLk∗ ‖2] + γ2kM2 +Ak.
By summing up k from 0 to n,
2 n∑ k=0 γkE[g(θk;Lk)− g(θLk∗ ;Lk)] ≤ −E[‖θn+1 − θ Ln+1 ∗ ‖2 − ‖θ0 − θ0∗‖2] +M2 n∑ k=0 γ2k + n∑ k=0 Ak
≤ D2 +M2 n∑ k=0 γ2k + n∑ k=0 Ak, (27)
where D = ‖θ0 − θ0∗‖. Dividing 2 ∑n k=0 γk on both sides, we have 1∑n k=0 γk n∑ k=0 γkE[g(θk;Lk)− g(θLk∗ ;Lk)] ≤ D2 +M2 ∑n k=0 γ 2 k + ∑n k=0Ak 2 ∑n k=0 γk ≤ D2 +M2 ∑n k=0 γ 2 k +A 2 ∑n k=0 γk .
According to the convexity of g(·; ·) and Jensen’s inequality Jensen (1906),
1∑n k=0 γk n∑ k=0 γkE[g(θk;Lk)] ≤ E[g(θ̄n; L̄n)], (28)
where θ̄n = ∑n
k=0 γkθk∑n k=0 γk
and L̄n = ∑n
k=0 γkLk∑n k=0 γk .
Then we have E[g(θ̄n; L̄n)]− ∑n k=0 γkg(θ Lk ∗ ;Lk)∑n
k=0 γk ≤ D2 +M2
∑n k=0 γ 2 k +A
2 ∑n k=0 γk
, (29)
We choose γk such that ∑n k=0 γk = ∞ and ∑n k=0 γ 2 k < ∞, for example, γk = 1k . Taking n to infinite, we have
0 ≤ lim n→∞ E[g(θ̄n, t̄n)]− lim n→∞
∑n k=0 γkg(θ
Lk ∗ , Lk)∑n
k=0 γk ≤ 0. (30)
Because Lk ↗ 0, g(·, ·) is continuous and ∑n k=0 γk =∞, then
lim n→∞
∑n k=0 γkg(θ
Lk ∗ , Lk)∑n
k=0 γk = g(θ0∗, 0), (31)
which implies that lim n→∞ E[g(θ̄n, L̄n)] = g(θ0∗, 0). (32)

heading: 4 NUMERICAL EXPERIMENTS

heading: 4.1 A NON-CONVEX EXAMPLE
text: We first use the following multi-dimensional highly non-convex example Chow et al. (2013) to illustrate the feasibility and efficiency of the HTA.
min f(x) = (π/n) { k sin2(πy1) +
n−1∑ i=1
(yi −A)2[1 + k sin2(πyi+1)] + (yn −A)2 } , (33)
where
x = (x1, · · · , xn) ∈ Rn, yi = 1 + (xi − 1)/4, i = 1, · · · , n and k = 10, A = 1. (34) The function has roughly 5n local minima in the region {|xi| < 10} and a unique global minimum located at xi = 1, i = 1, · · · , n. (See Figure 4.1 for an illustration of n = 2.) The idea of the homotopy setup of this problem is splitting the n dimensional optimization problem into two n/2 dimensional problems which can be optimized separately and are tracked to an optimal of solution of f(x) with respect to the homotopy parameter L. Then the objective function of the homotopy method is written as follows:
F (x, y;L) = f(x) + L‖(x1, · · · , xn2 )− (x 0 1, · · · , x0n2 )‖ 2 + g(y) + L‖(yn 2 +1 , · · · , yn)− (y0n 2 +1 , · · · , y0n)‖2
+ 1
L ‖x− y‖2,
(35) where x = (x1, · · · , xn) ∈ Rn, y = (y1, · · · , yn) ∈ Rn, (36) (x01, · · · , x0n 2 ) ∈ Rn2 , (y0n 2 +1 , · · · , y0n) ∈ R n 2 are fixed and chosen based on the initial guess. When L is large, F is equivalent to optimize
min xn
2 +1,··· ,xn g(x01, · · · , x0n2 , xn2 +1, · · · , xn) + miny1,··· ,yn 2 g(y1, · · · , yn2 , y 0 n 2 +1 , · · · , y0n), (37)
which can be solved separately. When L becomes small, the last term of (35) forces x and y to be equal. Then F (x, y;L) is equivalent to optimizing f(x) which is the original n-dimensional optimization problem. Therefore F (x, y;L) connects two optimization problems by gradually decreasing penalty factor L and builds a continuous path between the decoupled systems and the original optimization function f(x).
We compared the traditional optimization method (the quasi-Newton method) and the quasi-Newton with the homotopy setup. For n = 4, we chose the same initial value [5, 4, 5, 5] and set L from 0.9 to 0.1 for the quasi-newton with the homotopy setup. It takes 8 steps and 14 steps for the quasi-newton method and the quasi-newton with the homotopy setup to converge, respectively. However, the quasi-newton method does not find the global minimum while the HTA does. We also compared two algorithms with different initial values and count the times when converging to global minimums. Table 1 shows the results of n = 2, n = 4, n = 6 and n = 8.

heading: 4.2 APPLICATIONS ON VGG
text: Secondly, we applied the HTA to both CNN and fully connected neural networks on the CIFAR10 dataset (https://www.cs.toronto.edu/˜kriz/cifar.html) with the training and validation datasets from the CIFAR-10 website. In terms of the pre-processing step: for training dataset, we used random horizontal flip, random crop and normalization; for the validation dataset, we just used the same normalization. We implemented the HTA with stochastic gradient descent (SGD) method on VGG11, VGG13, VGG16 and VGG19 as our base models with batch normalization Simonyan & Zisserman (2014). In order to compare the SGD with the HTA method, we set all the regular parameters (e.g. learning rate, batch size, momentum, epochs, weight decay and etc) to be exactly the same, varied the homotopy parameter L from L = 0.01 to L = 0.005 with a stepsize ∆L = 0.001. For each L, we run 4 epochs with the SGD method. The code for this method can be found at https://github.com/Bill-research/homotopy which has been run on the DGX-1 (p100 version) with 8 p100 GPUs and Tesla structures. Figure 2 shows the comparison of validation loss between the HTA with fully connected neural networks and traditional method on the VGG13 model. The HTA has a less error rate (5.14%) than the traditional method (6.35%) by improving 11.68%. All the results for different models are shown in Table 2. It is clearly seen that accuracies of HTA for different models are better than that of the traditional method. For example, vgg11 with HTA gets 7.02% error rate while traditional method only gets 7.83%; The vgg16 with HTA has an error rate of 5.88% while traditional method only has 6.35%. We run the numerical experiments for 5 times for each model and the results are shown in Table 3. We also applied HTA to the only CNN part of two VGG models: vgg11 and vgg13. Figure 3 shows the validation accuracies of two models which has been improved by the HTA (91.61% (vgg11) and 93.09% (vgg13)) comparing with the traditional method (91.01% (vgg11) and 92.29% (vgg13)).

heading: 5 CONCLUSION
text: In this paper, we developed a homotopy training algorithm for solving the optimization problem arising from neural networks. This algorithm starts from decoupled low dimensional systems and gradually transforms to the coupled original system with complex structure. Then the complex neural networks can be trained by the HTA with a better accuracy. The convergence of the HTA for any given L is proved for the non-convex optimization. Then existence of solution path θ(L) is demonstrated theoretically for the convex case although that it exists numerically in the non-convex case. Several numerical examples have used to demonstrate the efficiency and feasibility of HTA. The application of HTA to VGG models on CIFAR-10 provides a better accuracy than the traditional method. Moreover, the HTA method would provide a new way to couple with the dropout techniques to improve the training accuracy which we will explore in the future.

Review summaries:
In this paper, the authors propose the Homotopy Training Algorithm (HTA) for neural network optimization problems. They claim that HTA starts with several simplified problems and tracks the solution to the original problem via a continuous homotopy path. They give the theoretical analysis and conduct experiments on the synthetic data and the CIFAR-10 dataset. 

The work proposes to learn neural networks using homotopy-based continuation method. The method divides the parameter space into two groups (extendable to multiple groups) and introduces a homotopy function which includes the original optimization problem as an extreme case.  By varying the homotopy parameter, one can construct a continuous path from a supposedly easier to solve optimization problem to the problem of interest. The authors prove convergence in the non-convex case, the existence of solution path in the convex case and demonstrate the effectiveness of the proposed method on synthetic and real datasets.