Andreas Griewank
email: griewank@math.hu-berlin.de

Ángel Rojas
email: angel.rojas@yachaytech.edu.ec

Bruce Watson

Eric Badouel

Eds. Oumar Niang

Generalized Abs-Linear Learning by Mixed Binary Quadratic Optimization

Keywords: Abs-Normal/Linear Form, Successive Piecewise Linearization, Mixed Binary Linear/Quadratic Optimization, Proximal Term. MOTS-CLÉS : Forme Abs-Normale / Linéaire, Linéarisation successive par morceaux, Optimisation binaire mixte linéaire / quadratique, Terme proximal

We consider predictor functions in generalized abs-linear form, which generalize neural nets with hinge activation. To train them with respect to a given data set of feature-label pairs, one has to minimize the average loss, which is a multi-piecewise linear or quadratic function of the weights, i.e. coefficients of the abs-linear form. We suggest to attack this nonsmooth, global optimization problem via successive piecewise linearization, which allows the application of mixed binary convex quadratic optimization codes amongst other methods. These solve the sequence of abs-linear model problems with a proximal term. Preliminary experiments on a simple regression problem verify the validity of the approach but require a large number of Simplex pivots by the solver Gurobi.

RÉSUMÉ.

Nous considérons les fonctions prédictives sous une forme abs-linéaire généralisée, qui généralisent les réseaux neuronaux avec activation de la charnière. Pour les entraîner par rapport à un ensemble de données donné de paires étiquette-caractéristique, il faut minimiser la perte moyenne, qui est une fonction linéaire ou quadratique multi-morceaux des poids, c'est-à-dire des coefficients de la forme abs-linéaire. Nous suggérons d'attaquer ce problème d'optimisation globale non lisse via une linéarisation successive par morceaux, qui permet l'application de codes d'optimisation quadratiques convexes binaires mixtes entre autres méthodes. Ceux-ci résolvent la séquence de problèmes du modèle abs-linéaire avec un terme proximal. Des expériences préliminaires sur un problème de régression simple vérifient la validité de l'approche mais nécessitent un grand nombre de pivots Simplex par le solveur Gurobi

Introduction and Notation

Neural nets with hinge activation have been proven theoretically [START_REF] Yarotsky | Error bounds for approximations with deep ReLU networks[END_REF] and experimentally [START_REF] Bölcskei | Optimal Approximation with Sparsely Connected Deep Neural Networks[END_REF] to produce prediction functions f (w; x) that are able to represent a wide variety of relations in machine learning. These predictor models are piecewise linear [START_REF] Scholtes | Introduction to Piecewise Differentiable Functions[END_REF] with respect to the feature vector x and multi-piecewise-linear (see explanation below) with respect to the weight vector w, which consists of various transformation matrices and inhomogeneous shifts. It is well known that every piecewise linear vector function from x ∈ R n to y ∈ R m can be expressed in an abs-linear form

y = f (w; x) ≡ b + Jx + N z + Y|z| s.t. z = F (w; x) = c + Zx + M z + L|z| (1)
where z ∈ R s is a vector of switching variables and the various coefficient vectors and matrices can be combined to the weight vector w ≡ (c, Z, M, L, b, J, N, Y ) ∈ R (s,s×n,s×s,s×s,m,m×n,m×s,m×s) R s with s = (s+m)(1+n+2s). To make sure that f (w; x) can be unambiguously evaluated as a piecewise linear continuous function of x we assume throughout that the square matrices M, L ∈ R s×s are strictly lower triangular.

In the case of multi-layer neural networks z represents nodal values and the matrices M and L are block diagonal, with M = L in the case of hinge activation. By induction on the s components z i of z one can easily see that they are piecewise linear functions z i (x) of x, which then also holds for the resulting vector

y = f (w; x) = b + Jx + N z(x) + Y |z(x)| ∈ R m .
However, it is important to note that the dependance of f (w; x) on the coefficients w is only multi-piecewise linear, i.e. piecewise linear with respect to each component of w when the others are kept constant. In other words for each Cartesian basis vector e j ∈ R s and fixed x the univariate function f (w + te j ; x) is piecewise linear. Therefore one can rather efficiently perform global coordinate searches even when the loss function on y is not piecewise linear but for example quadratic.

Separation, Scaling and Fixed Point Iteration

For the subsequent optimization approaches we normalize the abs-linear system by relegating the nonsmoothness to the second part of (1) and then scaling the z so that we can compute Lipschitz constants and an upper bound for z(x).

Separation of Nonsmoothness:

The smoothness condition Y = 0 can always be achieved by extending the originally given switching vector z to z = (z , |z| Y ) and correspondingly extending L, M and N to

L = L 0 Y 0 ∈ R (s+m)×(s+m) M = M 0 0 0 and Ñ = N I ∈ R m×(s+n) .
Naturally, we also have to pad the vector c ∈ R s and the matrix Z ∈ R s×n by m additional zero rows. Now the dimension of the coefficient vector w ≡ (c, Z, M, L, b, J, N ) has gone down by ms to (s + m)(1 + n) + s(m + 2s) after the incrementation of s of course. Thus a typical empirical risk function may look like

ϕ(w) ≡ 1 2 k k k=1 f (w, x k ) -ỹk 2 2 with f (w; x) = N z(w; x) (2)
where the (feature, label) pairs (x k , ỹk ) ∈ R n×m for k = 1, 2, . . . k form a suitable training set. Here N is typically a projection onto the last components of z.

Scaling for Contraction:

Secondly we can scale z to z = Dz with D = diag(d) a matrix of positive entries 0 < d i for i = 1 . . . s. Then the state equation can be rewritten as

Dz ≡ z = c + Zx + M z + L|z| ≡ Dc + DZx + (DM D -1 )Dz + (DLD -1 )|Dz| .
Hence we see that the strictly lower triangular matrices M and L undergo a positive diagonal similarity transformation. Starting from d i = 1 the d i can be chosen for i = 2 . . . s such that the rows of M and L have 1 norms less than or equal to any given bound ε 2 with ε < 1. This implies for their ∞ matrix norms

ε ≥ M ∞ + L ∞ ≥ | M | + | L| ∞ = M , L ∞
Obviously, it is sufficient that the sums of the 1 norms of corresponding rows in M and L, which equal to the rows obtained by summing or concatenating |M | and |L| are less than or equal to ε. Note that if scaling by D is performed, then the matrix N in the response must be transformed to ND -1 . Assuming without loss of generality that the norm bounds already hold for the original M and L we find that the map

G(z) ≡ G (M,L) (z) ≡ z -H(z) with H(z) = M z + L|z| (3)
is a perturbation of the identity by the map H(z). The latter has the componentwise Lipschitz property

|H(z) -H(z)| = |M (z -z) + L(|z| -|z|)| ≤ |M ||z -z| + |L|||z| -|z|| ≤ |M ||z -z| + |L||z -z| = (|M | + |L|)|z -z| ∈ R s (4)
where we have used that by the inverse triangle inequality ||z| -|z|| ≤ |z -z| componentwise. This implies for the infinity norms of matrices and vectors

H(z) -H(z) ∞ ≤ M, L ∞ z -z ∞ ≤ ε z -z ∞ ∈ R
so that we have a contraction provided ε < 1. Thus it follows as a consequence of the Banach fixed point theorem that there is a unique point z * = c + Zx + H(z * ), and G has an inverse G -1 with the Lipschitz constant 1/(1 -ε) < ∞. Hence we obtain for the exact solution z * the bound

z * = G -1 (M,L) (c + Zx) =⇒ z * ∞ ≤ G -1 (M,L) (c) ∞ + Z ∞ x ∞ (1 -M, L ∞ ) . ( 5
)

Fixed Point Iteration:

Moreover starting from z 0 we obtain for the iteration

z k+1 = c + Zx + M z k + L|z k | = c + Zx + H(z k ) (6)
from ( 4) the componentwise bound

|z k -z * | ≤ (|M | + |L|)|z k-1 -z * | ≤ (|M | + |L|) k |z 0 -z * | . (7)
In fact we do not only have a contraction but can establish finite convergence as follows.

Since |M | + |L| is strictly lower triangular we know that there is a minimal integer ν such that

(|M | + |L|) ν = 0 with ν ≤ s and ν = 0 ⇐⇒ M = 0 = L .
We call ν the switching depth since it counts how often the nonmooth elemental abs() can be super imposed on each other. In terms of the evaluation graph [START_REF] Griewank | Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation[END_REF] it is the longest directed path along abs() nodes. If ν = 0 the prediction function is smooth and thus linear, and if ν = 1 we say that the problem is simply switched, which is typical for KKT like systems. In the case of neural networks ν is simply the number of intermediate layers.

Because of the bound [START_REF] Bölcskei | Optimal Approximation with Sparsely Connected Deep Neural Networks[END_REF], the iteration (6) must reach z k = z * after at most ν steps from any z 0 . This applies irrespective of the scaling. However, in general only the norm D(z k -z * ) ∞ , and not the norm z k -z * ∞ , will decline monotonically. Of course from a good warm start we may hope to converge to a satisfactory accuracy in fewer than ν steps. We also emphasize that a matrix vector product can be much better computed in parallel rather than forward and backward substitutions on triangular matrices.

Mixed Binary Linear/Quadratic Optimization

For fixed w the constraints in (1) are essentially linear, except for the absolute values. Using general nonsmooth optimization methods in this situation would of course disregard a lot of structure. We will reformulate the problem as a mixed binary linear optimization problem, for which standard solvers are available.

Formulation w.r.t. x ∈ R n :

Assuming without loss of generality that Y = 0 as justified in Section 2 we can concentrate on the triangular state equation

z = c + Zx + M z + L|z| = c + Zx + M z + La with a = |z| ∈ R s . (8)
In order to eliminate the highly nonlinear constraint a = |z| we use a signature matrix Σ = diag(σ) whose elements σ i ∈ {-1, 1} for i = 1 . . . s are binary variables. Then the components of the vector a = Σz = (σ i z i ) i=1...s are bilinear in the binary variable σ i and the real variables z i . It is well known [START_REF] Gupte | Solving Mixed Integer Bilinear Problems Using MILP Formulations[END_REF] that the constraints a = Σz ≥ 0 can be linearized to the system of inequalties

-a ≤ z ≤ a and a + γ(σ -e) ≤ z ≤ -a + γ(σ + e) , (9)
where e ∈ R s is the vector of ones and γ ∈ R an upper bound on the possible norm values z ∞ . It can be computed by maximizing the right hand side of ( 5) over a range of x that are likely to occur.

From the property γ ≥ z ∞ it can be easily derived that these four vector inequalities are equivalent to the nonlinear constraint a = |z|. The key property of the inequalities ( 9) is that together with the equation z = c + Zx + M z + La it can be entered into Gurobi and other modern mixed integer optimization solvers to characterize a feasible set with respect to the real variable vector (x, z, a) ∈ R n+2s and the binary variable vector σ ∈ {-1, 1} s .

The corresponding objective must be linear or convex quadratic with respect to the real variables so that for example

f (w; x) = b + Jx + N z(x) + 1 2 (x -x) Q(x -x) ∈ R (10)
with some positive semi-definite matrix Q ∈ R n×n and a reference point x ∈ R n . We can also include quadratic terms with respect to z, which we will do later in the learning context. Frequently we may simply have an Euclidean proximal term in that Q = qI for some scalar 0 ≤ q ∈ R. Given any parameter q > 0 the objective will be bounded below and a single call to Gurobi will yield a global minimizer.

Proximal Term via Quadratic Overestimation

One possible origin of such a proximal term is that we have originally a nonlinear state equation

z = F (x, z, |z|) with F ∈ C 2 (R (n+2s) , R s ) . (11)
It can be abs-linearized at a reference point x with z = z(x) to

z = z + Z(x -x) + M (z -z) + L(|z| -|z|) (12)
where the matrices

Z = ∂ ∂x F (x, z, |z|) ∈ R (s×n) , M = ∂ ∂z F (x, z, |z|) ∈ R (s×s) L = ∂ ∂|z| F (x, z, |z|)
are all evaluated at the point (x, z, |z|). Again the square matrices M , L of order s must be structurally strictly lower triangular so that z(x) and z(x) can be evaluated unambiguously by ( 11) and (12), respectively. As shown in [START_REF] Griewank | On Stable Piecewise Linearization and Generalized Algorithmic Differentiation[END_REF] the error in this abs-linearization can be bounded by z(x) -z(x) ∞ ≤ q 2 x -x 2 for some q > 0 .

If the original scalar response was just f (w; x) = b+Jx+N z(x) we get the upper bound

f (w; x) ≤ b + Jx + N z(x) + N ∞ z(x) -z(x) ≤ f (w; x) ≡ b + Jx + N z(x) + q 2 x -x 2 with q = q N ∞ . ( 13
)
Hence we have exactly an objective of the proxlinear form (10) with z(x) defined by the abs-linear triangular system (12). It can of course be rewritten as a mixed binary linear system using the reformulation ( 9) and then solved by Gurobi.

While upper bounds for q = q N ∞ can be derived from the evaluation procedure of F at (x, z, |z|) it may also be updated iteratively like the size parameters in trust region or quadratic overestimation methods. As usual being conservative, i.e. making q large enforces convergence but may slow down the actual computation.

Outer Loop:

So the full method consists of forming the constraints (12) at the current iterate x, replacing |z| with an a satisfying ( 9) and minimizing the objective (13) by Mixed Binary Quadratic Optimization. The global minimizer x can then serve as the new reference point x provided we have for the real objective b

+ J x + N z(x) < b + Jx + N z(x) .
Otherwise the estimate for q was too small and must be increased significantly without a change in x. As shown in [START_REF] Fiege | An Algorithm for Nonsmooth Optimization by Successive Piecewise Linearization[END_REF] for a method that only computes local minimizers of the successive abs-linearizations, all clusterpoints x * of the sequence generated in this way are first order minimal, i.e. minimizers of the abs-linearization at x * . Unfortunately, even globally solving each local abs-linear problem does in general not guarantee that the cluster point is a global minimizer of the underlying problem.

Piecewise Linearization w.r.t. w:

So far we have assumed that the coefficients w = (c, Z, M, L) are constant throughout the minimization. For abs-linear learning we cannot make this assumption and rather have to consider w as variable and x as constant. Now the vector constraint z = c + Zx + M z + La contains the terms M z and La which are bilinear in the components of the real variables (z, a) and w. That means unfortunately that they cannot be linearized exactly.

To overcome this difficulty we perform a piecewise linearization of the original state equation z = c + Zx + M z + L|z| at a reference point (c, Z, M , L, z) with x fixed and z the corresponding z value such that z = c + Zx + Mz + L|z|. Thus we obtain with G as defined in (3) immediately the relation

G(z) = G (M,L) (z) = c + Zx -∆Mz -∆L|z|
where naturally ∆M = M -M and ∆L = L-L. Assuming that M and L are scaled as described above and also M, L ∞ ≤ ε < 1 we find that the inverse operator

G -1 = G -1 (M,L) still has the Lipschitz constant 1/(1 -ε).
Then we obtain for ∆z = z -z the bound

∆z ∞ ≤ ∆c + ∆Zx + ∆Mz + ∆L|z| ∞ 1 -M, L ∞ ≤ ∆c ∞ + ∆Z ∞ x ∞ + ∆M, ∆L ∞ z ∞ 1 -M, L ∞ (14)
where naturally ∆c = c -c and ∆Z = Z -Z. Now we do want to find a piecewise linearization of the abs-normal function

z = F (w) = G -1 (M,L) (c + Zx)
i.e. an abs-linear function z = F (w) approximation in the variables w = (c, Z, M, L).

Proposition. The vector function

z = F (w) ≡ G -1 ( M , L) (c + Zx + ∆Mz + ∆L|z|)
is abs-linear and can be evaluated by solving the triangular system

z = (c + Zx + ∆Mz + ∆L|z|) + M z + L|z| . ( 15
)
Moreover, provided M , L ∞ + ∆M, ∆L ∞ ≤ ε < 1 (16)
the discrepancy between F (w) and its abs-linearization F (w) is bounded by

F (w) -F (w) ∞ ≤ z ∞ (1 -ε) 2 ∆M, ∆L 2 ∞ ≤ s z ∞ (1 -ε) 2 ∆M, ∆L 2 F .
The proof of this result can be found in the appendix Section 6.

The overall Mixed Binary Quadratic Problem:

Using (z+ z) N N (z-z) = N z 2 -N z 2
we obtain for a single loss term in (2)

1 2 N z-ỹ 2 2 -1 2 N z-ỹ 2 2 = 1 2 (N (z + z)-2ỹ) N (z-z) ≤ Nz-ỹ 2 N (z-z) 2 +O( ∆w 3 2 ) ≤ Nz-ỹ 2 N 1 s z ∞ (1 -ε) 2 ∆M, ∆L 2 F +O( ∆w 3 2 ) . (17)
For each one of the sample points (x k , ỹk ) we get a different z = z(w; x k ) which will be denoted simply by z k (w). Then the resulting mixed binary quadratic problem has the convex quadratic objective

ϕ(w) ≡ 1 2 k k k=1 N z k (w) -ỹk 2 2 + s N 1 ∆w 2 2 k(1 -ε) 2 k k=1 zk ∞ Nz k -ỹk 2 . ( 18
)
So the first term is the approximation to the original quadratic loss risk via the piecewise linearization of the state equation and the second term is a proximal quadratic that tries to bound the resulting error in the objective.

For each data index k = 1 . . . k we get with a k ∈ R s and σ k ∈ {-1, +1} s the bilinear constraint set

z k = c + Zx + ∆Mz k + ∆L|z k | + M z k + La k (19) -a k ≤ z k ≤ a k and a k + γ(σ k -e) ≤ z k ≤ -a k + γ(σ k + e) .
Hence we see that we have the same binary variables as in the formulation w.r.t. x and the only thing that changes compared to ( 8) is that we have the new linear terms ∆Mz and ∆L|z|. The bound γ can of course be adjusted individually to γ k for each k. The minimizing of (18) subject to the linear constraints ( 19) and ( 9) with the σ k binary can be delegated to solvers like Gurobi. The global minimizer ŵ of the model problem can then serve as the reference point ẘ for the next abs-linearization, where we have to recompute the proximal term in (18).

Preliminary Numerical Experiment and Conclusion

Given the limitations of space and time we have not been able to conduct and report numerical results that are in any way conclusive. We simply minimized the averaged loss defined in (2) for the Griewank function [START_REF] Griewank | Generalized descent for global optimization[END_REF] ỹ(x) = 1 + 1 4000

  over 50 training points and 8 testing points x ∈ R 2 chosen uniformly at random in the cube [-8, 8] d . The results reported are for the model sizes s = 3, 4, 5 over 5 successive linearizations. This results in 15 mixed integer quadratic optimization problems specified in (2) and (19). The following graph shows the true objective function values achieved by successive linearizations.As one can see for the two smaller prediction models the minimal objective function is essentially already reached at the second linearization whereas the model with s = 5 switching variables can bring the objective further down. The number of weights, Gurobi variables and Simplex iterations are listed in the following table. The bad news is that the latter appears to be growing rapidly with respect to the model size. Nevertheless the exact solution can serve as a reference point for other methods and hopefully Gurobi can be accelerated by exploiting more of the special structure.

				d k=1	x 2 k -	d k=1	cos	x k √ k
	0.16
	Empirical Risk							Model size
								s=3
	0.15							s=4
								s=5
	0.14
	0.13
	0.12
	0.11
								# of models
	1	2			3			4	5
	s #w #variables	pl 1		pl 2		pl 3	pl 4	pl 5
	3 21	471	303810	353703	1716277	581060	681025
	4 31	631 1129639	263007	1015447	1339147	1068608
	5 43	793 1153345 22793377 22895320 21241422 16513124

Proceedings of CARI 2020

Appendix

Proof(Proposition). According to [START_REF] Fiege | An Algorithm for Nonsmooth Optimization by Successive Piecewise Linearization[END_REF] the exact solution z = F (w) satisfies the triangular system z -M z -L|z| = c + Zx + ∆M z + ∆L|z| and thus

The difference between the two right hand sides of z and z is ∆M (z -z) + ∆L(|z| -|z|) so that by the Lipschitz continuity of G -1

we can collect the terms in z -z ∞ on the left hand side and then divide both sides by its factor yielding

The penultimate inequality follows since again by Lipschitz continuity of G -1 ( M , L) between the arguments c + Zx and c

For the final bound we use that the the ∞ norm of a (2s × s) matrix is the maximal 1 norm of any one of its s rows, which is the bounded by √ 2s times its 2 vector norm, which in turn is bounded by the Frobenius norm of the full matrix. This bound is sharp when ∆M, ∆L has one nontrivial row of constants.