Year: 2019

Abstract:
Preconditioning an minimization algorithm improve its convergence and can lead to a minimizer in one iteration in some extreme cases. There is currently no analytical way for finding a suitable preconditioner. We present a general methodology for learning the preconditioner and show that it can lead to dramatic speed-ups over standard optimization techniques.

heading: 1 INTRODUCTION
text: Many problems arising in applied mathematics can be formulated as the minimization of a convex function f : Rd → (−∞,+∞]
min x∈Rd f(x).
The resolution of an optimization problem is usually tackled using an optimization algorithm that produces a sequence of iterates converging to some minimizer of f (Nesterov (2018)). The gradient descent algorithm is a standard optimization algorithm that converges linearly (i.e exponentially fast) if f is regular enough. Preconditioned methods (Nemirovsky & Yudin (1983); Lu et al. (2016); Maddison et al. (2019)) are powerful optimization algorithms that converges linearly under weaker assumptions that gradient descent. The performance of these methods relies heavily on the preconditioning – the task of choosing the hyperparameter called the preconditioner. In the case of the Dual space Preconditioned Gradient Descent (DPGD) of Maddison et al. (2019), an optimal preconditioning can lead to convergence in one iteration.
The preconditioner is a function that has to be selected properly w.r.t. f to obtain the desired linear convergence. Although Maddison et al. (2019) gives some hints to precondition DPGD, there is currently no analytical way for finding a suitable preconditioner.
In this paper, we propose to learn the preconditioner of DPGD using a neural network. We make the following contribution:
• We propose a supervised learning setting to learn the preconditioner of an optimization algorithm (DPGD, Maddison et al. (2019))
• We present a general methodology that allows to effectively learn the preconditioner while avoiding issues related to this task
• We implement this methodology in dimension one and 50 and show that this can lead to dramatic speed-ups.
The remainder is organized as follows. The next section provides background knowledge on the DPGD algorithm. Then, in section 3 we present our supervised learning setting. Finally, we apply this methodology in section 4. Additional developments are provided in the appendix as well as postponed proofs.

heading: 2 BACKGROUND ON DUAL PRECONDITIONING
text: The main message of this section is that a good preconditioner is a preconditioner that is "close" to ∇f?, where f? is the Fenchel conjugate of f .

heading: 2.1 LEGENDRE FUNCTIONS
text: Before introducing DPGD algorithm, we need to define the class of functions f that we aim to minimize. These functions are called Legendre convex functions and are studied in (Rockafellar, 1970, Section 26). We recall their definition and some of their properties. Given a function f : Rd → (−∞,∞], the domain of f is the set {x ∈ Rd, f(x) < ∞}. We denote D(f) the interior of the domain of f .
Definition 1. The function f is Legendre if D(f) is not empty, and if
• the function f is lower semicontinuous on Rd, and is differentiable and strictly convex on D(f),
• the gradient∇f is coercive on D(f): For every sequence (xi)i∈N ∈ D(f) converges to a point on the boundary of the domain of f , ‖∇f(xi)‖ → ∞.
Given a Legendre function f , one can define its Fenchel conjugate f?(y) = supx∈Rd〈x, y〉 − f(x). Proposition 1 (Rockafellar (1970)). The function f is Legendre if and only if its Fenchel conjugate f? is Legendre, in which case ∇f is one-to-one between D(f) and D(f?). Moreover, ∇f is continuous in both directions, for every x ∈ D(f), ∇f?(∇f(x)) = x and for every y ∈ D(f?), ∇f(∇f?(y)) = y. Proposition 2 (Rockafellar (1970)). If f is Legendre, then it either has no minimizer, or one unique global minimizer x? ∈ D(f). If f admits a minimizer, then ∇f(x?) = 0, 0 ∈ D(f?) and x? = ∇f?(0).

heading: 2.2 DPGD ALGORITHM (MADDISON ET AL. (2019))
text: Consider a convex Legendre function f : Rd → (−∞,+∞]. To minimize f , the DPGD algorithm is written x+ = x− γ [∇p(∇f(x))−∇p(0)] , (1) where γ > 0 and p : Rd → (−∞,+∞] are parameters. DPGD is well defined over the set of Legendre convex functions f , see Maddison et al. (2019). The map ∇p is called a non-linear preconditioner of ∇f . The preconditioning is the design of∇p.

heading: 2.3 CONVERGENCE THEORY OF DPGD
text: In the sequel, we consider a Legendre function f admitting a minimizer x?. The preconditioning can dramatically affect the convergence of DPGD. Indeed, if ∇p = ∇f?, then DPGD algorithm converges in one iteration. To see this, note that in this case
x+ = x− [x−∇p(0)] = ∇p(0) = ∇f?(0),
if γ = 1 and recall that∇f?(0) is the minimizer of f (Proposition 2). Of course, this preconditioning is unrealistic since computing∇f?(0) is as hard as minimizing f . It is more realistic to aim to find a preconditioner p that satisfies the following conditions.
Assumption 1 (Dual relative smoothness). There exists L ≥ 0 such that
〈∇p(x)−∇p(y), x− y〉 ≤ L〈∇f?(x)−∇f?(y), x− y〉, ∀x, y ∈ D(p). (2)
Assumption 2 (Dual relative strong convexity). There exists µ > 0 such that
〈∇p(x)−∇p(y), x− y〉 ≥ µ〈∇f∗(x)−∇f∗(y), x− y〉, ∀x, y ∈ D(p). (3)
If p and f? are twice differentiable, these assumptions can be written as
µ∇2f? (3) ∇2p (2) L∇2f?,
for the ordering of nonnegative matrices. These assumptions mean that∇p must be chosen close to∇f? in some sense. Moreover, they are sufficient for the linear convergence of DPGD.
Theorem 1 (Maddison et al. (2019), Informal). If assumption 1 holds and if p is Legendre, then the sequence of iterates of DPGD converges to x? = argmin f . If assumption 2 holds, then the convergence is linear.
In comparison, gradient descent (and many other first order methods, Nesterov (2018)) converges for smooth convex functions f (i.e. functions f such that ∇f is Lipschitz continuous). Moreover, the convergence is linear if f is strongly convex. Many functions satisfy assumption 1 (and assumption 2) for some preconditioner p, without being smooth (or strongly convex), see Maddison et al. (2019).
If f is smooth and strongly convex, the convergence rate of gradient descent depends on the global condition number κ of f (defined as the Lipschitz constant of ∇f divided by the strong convexity parameter). If κ is high, then the linear convergence of gradient descent is slow and f is said ill-conditioned. The function f is said well-conditioned else. When f is ill-conditioned, it might be well-conditioned locally (i.e., on small subsets of Rd) but the rate of convergence of gradient descent only takes into account the global, high condition number. This is not the case for preconditioned algorithms ( Li & Malik (2016); Maddison et al. (2019)).

heading: 3 SUPERVISED LEARNING SETTING
text: In this section, we describe our methodology to learn∇f?.

heading: 3.1 FORMULATION
text: Our idea is to sample points xi ∈ D(f) according to some distribution and create the dataset {(∇f(xi), xi), i ∈ {1, . . . , n}} where the ∇f(xi) represent the features and xi the labels. The feature space is D(f?) and the label space is D(f) (Proposition 1). Since∇f? maps∇f(xi) to xi (Proposition 1), ∇f? can be seen as the solution of this supervised learning problem. We choose a neural net to modelize∇f?. To solve our supervised learning problem, we design a machine learning algorithm. The goal of the training of the machine learning algorithm is to minimize the theoretical risk Lµ(θ) Lµ(θ) = EX∼µ (`(Model∇f∗(∇f(X), θ), X)) , (4) where µ is the distribution of the labels X , θ represents the parameters of the neural network, Model∇f?(∇f(X), θ) is the output of the neural network with input∇f(X) and parameter θ, and
`(Model∇f∗(∇f(X), θ), X)
is the loss associated to one sample (∇f(X), X). To fully specify the machine learning algorithm, we need to fully define Lµ and the way we approximate it. We proceed by
• Choosing an architecture to specify the function Model∇f∗ (this is problem specific, see sections 3.2, 4.2)
• Giving a value for µ (section 3.3) • Providing a way to approximate the expectation in (4) by a finite sum (section 3.4) • Choosing a loss ` (this is problem specific, see section 4.2).
Once this is done, the algorithm is trained using Stochastic Gradient Descent (SGD) (Bottou et al. (2018)).

heading: 3.2 GENERAL SETUP
text: In this section we give a general description of the map Model∇f∗ . According to recommendations in Bengio (2012), we first standardize the features∇f(xi) and the labels xi. In other words, we apply a diffeomorphism (i.e, a one-to-one map which is continuously differentiable in both directions) G : D(f?)→ [−0.5, 0.5]d to the features∇f(xi) and a diffeomorphism H−1 : D(f)→ [−0.5, 0.5]d to the labels xi. Then a neural network is used to learn a predictor from the standardized features to the standardized labels, see figure 1. One can obtain a label from a standardized label by applying the
inverse map H of H−1, which can be easily computed (in practice H−1 is essentially an addition and a multiplication by a positive number), see section 4.
3.3 CHOICE OF µ
We found that the choice of the distribution µ that parametrized Lµ is critical for the performance of the machine learning algorithm. The reason is the following. Recall that we want to learn ∇f? to precondition DPGD algorithm whose goal is to minimize nonsmooth and non-strongly convex (or ill-conditioned) objectives f , see section 2. Let us imagine that we chose an uniform distribution for µ. Then, the situation in dimension one is represented by the figure 2.
The image distribution ν of µ by ∇f over-represents some areas of the feature space while other areas are under-represented. When the distribution of the features is as degenerated as ν, it is obvious that machine learning algorithms cannot generalize well (see Bengio (2012); Mesnil et al. (2011)). Therefore, we have to choose µ carefully. The classical recommendation is to ensure that the input of the neural network (i.e the standardized features) are uniformly distributed (Mesnil et al. (2011)). In other words, we have to find µ such that, if xi ∼ µ, then G ◦ ∇f(xi) is uniformly distributed over [−0.5, 0.5]d. This can be done by using the change of variable formula.
Proposition 3. Assume that ∇f and ∇f? are differentiable. If µ is the distribution with density proportional to
|det(JG◦∇f (x))|
with respect to Lebesgue measure (where J denotes the Jacobian matrix), then the distribution of G ◦ ∇f(X) where X ∼ µ is uniform over [−0.5, 0.5]d.
In the sequel, µ is the distribution whose density is given by Proposition 3. Inspecting figure 2, this seems reasonable since we expect X ∼ µ to have a high density when∇f ’s slope is large and a low density when∇f ’s slope is small. Remark 1. It is worth noting that the evaluation of µ doesn’t require any inversion as in the inverse transform sampling technique. Besides, second order optimization algorithms usually require (an approximation of) the inverse of the Hessian matrix, which we don’t need.

heading: 3.4 ESTIMATION OF THE LOSS
text: At this step, µ is chosen and Lµ has to be minimized by SGD. As usual, Lµ cannot be computed in closed form. The classical idea is to approximate Lµ by an empirical mean involving samples from the distribution µ. Sampling from µ turns out to be difficult in our case. Indeed, the use of rejection-based MCMC methods lead to a high rejection rate. Since we can evaluate the density of µ, we suggest to use importance sampling instead. Proposition 4. If xi ∼ π where the density of π is positive, then
E
( 1
n n∑ i=1 `(Model∇f∗(∇f(xi), θ), xi) µ(xi) π(xi)
) = Lµ(θ),
where µ(x) (resp. π(x)) denotes the density of µ (resp. π) with respect to Lebesgue measure.
We call the map θ 7→ 1n ∑n i=1 `(Model∇f∗(∇f(xi), θ), xi) µ(xi) π(xi)
the empirical risk. Proposition 4 states that if xi ∼ π, then the empirical risk is an unbiased estimator of the theoretical risk. Since this estimator will be minimized using SGD, we only need to know it up to a constant factor. Therefore we can get rid of the 1/n factor and of constant factors defining the densities of µ and π. The distribution π is called the proposal distribution.

heading: 3.5 RELATED WORKS
text: Learning an optimization algorithm has already been proposed in the literature. In Li & Malik (2016); Andrychowicz et al. (2016) the authors propose to learn a better optimization algorithm by observing its execution over a feature space of objective functions. This outperforms existing hand-engineered algorithms in terms of convergence speed, but require to know the performance of some algorithms over some objectives. Using optimization literature knowledge, our approach allows to reduce the problem to a supervised learning setting in Rd, that we can solve with a neural net. In its philosophy, our approach can also be related to deep learning methods for inverse problems used in signal processing, see e.g Lucas et al. (2018).

heading: 4 NUMERICAL EXPERIMENT
text: In this section, we run DPGD algorithm after learning ∇f?. We first consider a one dimensional toy problem with a ground truth involving power functions as in Maddison et al. (2019). Then ill-conditioned logistic losses are considered. These functions cannot be minimized (efficiently) by gradient descent. We call LDPGD our approach (for Learned DPGD).

heading: 4.1 POWER FUNCTIONS AND LOGISTIC LOSS
text: Two class of objective functions are considered. Power functions are defined by
f : x 7→ 1 a ‖Ax− c‖a , where  a ∈ (2,+∞) c ∈ Rd A is a matrix.
(5)
These functions are nonsmooth and non-strongly convex, therefore gradient descent is not guaranteed to minimize them. Power function minimization is the main example of Maddison et al. (2019).
We also consider the logistic loss
f : x 7→ 1 n n∑ i=1 log (1 + exp (−bi 〈ai, x〉)) + r 2 ‖x‖2, (6)
where r > 0 and for every i, bi ∈ {−1, 1}n and ai ∈ Rd. These functions are strongly-convex and smooth. It is easy to control the condition number of f :
Lemma 1. Consider the logistic loss f and the matrix A whose lines are the ai, and denote λmax the largest eigenvalue of the real symmetric matrix ATA. Denote L the smoothness parameter of f 1 and κ the condition number of f . Then,
L ≤ L′ = λmax 4n + r, and κ ≤ λmax 4nr + 1.

heading: 4.2 EXPERIMENTAL SETUP
text: We use a three layers neural network, with leaky ReLU activations and respectively 256 and 128 neurons in the two hidden layers. The square loss is used `(y, x) = 12‖y − x‖
2. The proposal distribution π for the importance sampling is taken gaussian. In the case of power functions, the features are rescaled using the map
log-rescaling : x 7→ {
log(1 + x) if x > 0 − log(1− x) else.
Justification for this choice is provided in the appendix (section D). Then the rescaled features and the labels are normalized to have zero mean and unit variance. The general architecture is summarized in figure 3.
The two first blocks correspond to the map G and the last block to the map H of section 3.2.
The model is trained with a dataset of 1000 samples.

heading: 4.3 RESULTS FOR POWER FUNCTIONS
text: Consider f a power function (5) where d = 1, a = 50, c = 50 and A is randomly chosen. We train the model during 100 epochs.
The next figure represents the performance of the machine learning algorithm we developped to learn ∇f∗ (note the scale of the inputs of∇f∗, see section D).
1i.e, the Lipschitz constant of ∇f
We see that∇f? is accurately learned the machine learning algorithm. This approximation of∇f? is then used as a preconditioner of DPGD. The value of the objective function and the value of the iterates while running LDPGD is plotted in the next figure. The stepsize is set to 1 during 10 steps and LDPGD is initialized at x0 = 100.
The iterates of gradient descent always diverge (no matter the stepsize choice). On the contrary, LDPGD iterates quickly converge to 357.1 (with ∇f(357.1) = 3.8× 10−14) which is close to the actual minimizer (356.4 in this example).

heading: 4.4 RESULTS FOR LOGISTIC REGRESSION
text: Consider f a logistic loss (6) with d = 50, n = 1000, and condition number κ. We train the model during 100 epochs. and represent the evolution of the objective function and the iterates while running LDPGD (after learning∇f?) and gradient descent (GD). The algorithms are initialized at (50, 50, 50)T and we use a stepsize for which GD converges linearly2. This is done for several logistic losses with worsening condition numbers. To vary κ, we vary the parameter r that controls the strong convexity of the objective f (see Lemma 1). Intuitively, if r is high then f is well-conditioned and if r is low f is ill-conditioned. The left figure represents the iterates (each curve corresponds to the evolution of one coordinate of the iterates) and the right figure represents the objective fonction values while running both algorithms.
2GD converges linearly if the stepsize is smaller than 1/L, where L is the Lipschitz constant of ∇f . Therefore, it is enough to take a stepsize smaller than 1/L′ where L′ is defined in Lemma 1
Our careful preconditioning allows LDPGD to outperform GD especially when the objective is ill-conditioned, as predicted by the theory, see section 2.

heading: 1 Introduction 1

heading: 2 Background on dual preconditioning 1
text: 2.1 Legendre functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2.2 DPGD algorithm (Maddison et al. (2019)) . . . . . . . . . . . . . . . . . . . . . . 2
2.3 Convergence theory of DPGD . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

heading: 3 Supervised learning setting 3
text: 3.1 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
3.2 General setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
3.3 Choice of µ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
3.4 Estimation of the loss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
3.5 Related works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

heading: 4 Numerical experiment 5
text: 4.1 Power functions and logistic loss . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
4.2 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
4.3 Results for power functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
4.4 Results for logistic regression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

heading: A Proof of Proposition 3 11

heading: B Proof of Proposition 4 11

heading: C Proof of Lemma 1 11
text: D Log-rescaling 12

heading: A PROOF OF PROPOSITION 3
text: Recall that G is a diffeomorphism and that f is Legendre. Let U the uniform distribution over [−0.5, 0.5]d and denote q its density with respect to Lebesgue measure. Under the assumption of Proposition 3, ∇f : D(f) → D(f?) is a diffeomorphism. Therefore, the map M = G ◦ ∇f : D(f) → [−0.5, 0.5]d is also a diffeomorphism and its inverse map is M−1 = ∇f? ◦ G−1 : [−0.5, 0.5]d → D(f). Our problem is to find a distribution µ supported by D(f) such that M(X) = G ◦ ∇f(X) ∼ U if X ∼ µ. Such µ satisfies for every nonnegative measurable function φ,
∫ φ(x)dµ(x) = ∫ φ(M−1 ◦M(x))dµ(x)
= EX∼µ(φ(M−1 ◦M(X))) = EU∼U (φ(M−1(U)))
= ∫ φ(M−1(u))q(u)du
= ∫ D(f) φ(M−1(M(x)))q(M(x))|det(JM (x))|dx
= ∫ D(f) φ(x)q(M(x))|det(JM (x))|dx,
where JM denotes the Jacobian matrix of M and the penultimate equality follows from the change of variable formula (which is allowed because M is a diffeomorphism). Therefore, the only possible solution to our problem is the distribution µ with density proportional to |det(JM (x))| if M(x) ∈ [−0.5, 0.5]d (which is always satisfied by definition of M ) and 0 else. One can check that this density solves the problem.

heading: B PROOF OF PROPOSITION 4
text: We apply the importance sampling principle to the measurable nonnegative function (the parameter θ is fixed) φ(x) = `(Model∇f∗(∇f(x), θ), x). We have
EY∼µ(φ(Y )) = EX∼π ( φ(X) µ(X)
π(X)
) = E ( 1
n n∑ i=1 φ(xi) µ(xi) π(xi)
) ,
where xi ∼ π.

heading: C PROOF OF LEMMA 1
text: Denote L the smoothness constant of f (i.e the Lipschitz constant of its gradient) and λ its strong convexity parameter. For every x ∈ Rd,
∇f(x) = 1 n n∑ i=1 −biai 1 1 + exp (bi 〈ai, x〉)) + rx
∇2f(x) = 1 n n∑ i=1 b2i aia T i exp (bi 〈ai, x〉)) (1 + exp (bi 〈ai, x〉)))2 + rId,
where Id is the identity matrix. Note that for all a 6= 1, a(1+a)2 ≤ 1 4 , and that b 2 i ≤ 1. It follows that
∇2f A TA
4n + rId.
Therefore, L ≤ λmax4n +r. Since κ = L/λ, the result follows from the fact that r ≤ λ (see Equation 6).

heading: D LOG-RESCALING
text: For power functions,∇f(x) can range in [ −10100, 10100 ]d if x is sampled from a gaussian. Assume d = 1 that we use Doubles precision in the implementation. This means that all the values we manipulate are stored with a precision of approximately 15 decimal places.
If we naively standardize (∇f(xi))i∈[1,n] to have zero mean and unit variance, we rescale numbers in the range [−10100, 10100]d, to the range [−0.5, 0.5]d. This means that a precision of 10−15 in the standardized scale represents a precision of 1085 in the original scale.
In practice, this loss of precision means that inputs that are different before standardization, might be indistinguishable after the standardization. As a result, the neural network might receive many different outputs for a single input value.
As it tries to minimize the empirical risk, the model will tend to associate to input values an average of the corresponding outputs. As a consequence, the regressed model will be a piecewise constant function, and will be constant on subsets of diameter 1085.
This is especially problematic for inputs close to 0. Let us denoteN0 the neighborhood of 0 on which the regressed model is constant equal to Model∇f∗(0, θ). Any x ∈ D(f) for which∇f(x) ∈ N0 is a fixed point of the DPGD alogorithm:
x+ = x− 1 Lp [Model∇f∗(∇f(x), θ)−Model∇f∗(0, θ)]
= x− 1 Lp [0] = x
The consequence is very undesirable: the solution offered by the DPGD alogorithm might be x even though ∇f(x) = O(1085) is far from zero.
As mentioned above, the lack of precision due to a naive standardization of the features is especially problematic around 0. We therefore propose to preserve the scale around small values more than around large values. This can be achieved by applying the following transformation to the features before normalizing them (to have zero mean and unit variance):
log-rescaling : x 7→ {
log(1 + x) if x > 0 − log(1− x) else.

Review summaries:
This paper proposes an optimization method with the preconditioning in the framework of supervised learning. The ideal preconditioning is given by the Fenchel conjugate of the optimization function. This paper uses a supervised scenario to find the ideal preconditioning. The authors point out the importance of the sampling distribution then and propose a sampling scheme using the uniform distribution on the space of gradient descents. The samples are used to train neural networks that imitate the mapping of the Fenchel conjugate. The training network is incorporated into the Dual space Preconditioned Gradient Descent (DPGD). Some numerical experiments show the effectiveness of the proposed method comparing to the standard gradient descent method.