Next Article in Journal
Biologically Inspired Visual System Architecture for Object Recognition in Autonomous Systems
Next Article in Special Issue
On a Nonsmooth Gauss–Newton Algorithms for Solving Nonlinear Complementarity Problems
Previous Article in Journal
Image Edge Detector with Gabor Type Filters Using a Spiking Neural Network of Biologically Inspired Neurons
Previous Article in Special Issue
On the Use of Biased-Randomized Algorithms for Solving Non-Smooth Optimization Problems
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Polyhedral DC Decomposition and DCA Optimization of Piecewise Linear Functions

Institut für Mathematik, Humboldt-Universität zu Berlin, 10099 Berlin, Germany
*
Author to whom correspondence should be addressed.
Submission received: 28 May 2020 / Revised: 5 July 2020 / Accepted: 8 July 2020 / Published: 11 July 2020

Abstract

:
For piecewise linear functions f : R n R we show how their abs-linear representation can be extended to yield simultaneously their decomposition into a convex f ˇ and a concave part f ^ , including a pair of generalized gradients g ˇ R n g ^ . The latter satisfy strict chain rules and can be computed in the reverse mode of algorithmic differentiation, at a small multiple of the cost of evaluating f itself. It is shown how f ˇ and f ^ can be expressed as a single maximum and a single minimum of affine functions, respectively. The two subgradients g ˇ and g ^ are then used to drive DCA algorithms, where the (convex) inner problem can be solved in finitely many steps, e.g., by a Simplex variant or the true steepest descent method. Using a reflection technique to update the gradients of the concave part, one can ensure finite convergence to a local minimizer of f, provided the Linear Independence Kink Qualification holds. For piecewise smooth objectives the approach can be used as an inner method for successive piecewise linearization.

1. Introduction and Notation

There is a large class of functions f : R n R that are called DC because they can be represented as the difference of two convex functions, see for example [1,2]. This property can be exploited in various ways, especially for (hopefully global) optimization. We find it notationally and conceptually more convenient to express these functions as averages of a convex and a concave function such that
f ( x ) = 1 2 ( f ˇ ( x ) + f ^ ( x ) ) with f ˇ ( x ) convex and f ^ ( x ) concave .
Throughout we will annotate the convex part by superscript   ˇ and the concave part by superscript  ^ , which seems rather intuitive since they remind us of the absolute value function and its negative. Since we are mainly interested in piecewise linear functions we assume without much loss of generality that the functions f and the convex and concave components are well defined and finite on all of the Euclidean space R n . Allowing both components to be infinite outside their proper domain would obviously generate serious indeterminacies, i.e., NaN s in the numerical sense. As we will see later we can in fact ensure in our setting that pointwise
f ^ ( x ) f ( x ) f ˇ ( x ) for all x R n ,
which means that we actually obtain an inclusion in the sense of interval mathematics [3]. This is one of the attractions of the averaging notation. We will therefore also refer to f ^ and f ˇ as the concave and convex bounds of f.

Conditioning of the Decomposition

In parts of the literature the two convex functions f ˇ and f ^ are assumed to be nonnegative, which has some theoretical advantages. In particular, see, e.g., [4], one obtains for the square h = f 2 of a DC function f the decomposition
h = 1 4 ( f ˇ + f ^ ) 2 = 1 2 1 4 ( f ˇ 2 + f ^ 2 ) h ˇ + 1 4 [ ( f ˇ f ^ ) 2 ] h ^ .
The sign conditions of f ˇ and f ^ are necessary to ensure that the three squares on the right hand side are convex functions. Using the Apollonius identity f · h = 1 2 [ ( f + h ) 2 f 2 h 2 ] one may then deduce in a constructive way that not only sums but also products of DC functions inherit this property. In general, since the convex functions f ˇ and f ^ have both supporting hyperplanes one can at least theoretically always find positive coefficients α and β such that
f ˇ ( x ) + α + β x 2 0 f ^ ( x ) α β x 2 for x R n .
Then the average of these modified functions is still f and their respective convexity/concavity properties are maintained. In fact, this kind of proximal shift can be used to show that any twice Lipschitz continuously differentiable function is DC, which raises the suspicion that the property by itself does not provide all that much exploitable structure from a numerical point of view. We believe that for its use in practical algorithms one has to make sure or simply assume that the condition number
κ ( f ˇ , f ^ ) sup x R n | f ˇ ( x ) | + | f ^ ( x ) | | f ˇ ( x ) + f ^ ( x ) | [ 1 , ]
is not too large. Otherwise, there is the danger that the value of f is effectively lost in the rounding error of evaluating f ˇ + f ^ . For sufficiently large quadratic shifts of the nature specified above one has κ β . The danger of an excessive growth in κ seems akin to the successive widening in interval calculations and similarly stems also from the lack of strict arithmetic rules. For example doubling f and then subtracting it yields the successive decompositions
( 2 f ) f = ( f ˇ + f ^ ) 1 2 ( f ˇ + f ^ ) = ( f ˇ 1 2 f ^ ) + ( f ^ 1 2 f ˇ ) = 1 2 [ ( 2 f ˇ f ^ ) + ( 2 f ^ f ˇ ) ] .
If in Equation (3) by chance we had originally f ^ = 1 2 f ˇ > 0 so that f = 1 2 f ˇ with the condition number κ ( f ˇ , 0.5 f ˇ ) = 3 we would get after the doubling and subtraction the condition number κ ( 2.5 f ˇ , 2 f ˇ ) = 9 . So it is obviously important that the original algorithm avoids as much as possible calculations that are ill-conditioned in that they even just partly compensate each other.
Throughout the paper we assume that the functions in question are evaluated by a computational procedure that generates a sequence of intermediate scalars, which we denote generically by u , v and w. The last one of these scalar variables is the dependent, which is usually denoted by f. All of them are continuous functions u = u ( x ) of the vector x R n of independent variables. As customary in mathematics we will often use the same symbol to identify a function and its dependent variable. For the overall objective we will sometimes distinguish them and write y = f ( x ) . For most of the paper we assume that the intermediates are obtained from each other by affine operations or the absolute value function so that the resulting u ( x ) are all piecewise linear functions.
The paper is organized as follows. In the following Section 2 we develop rules for propagating the convex/concave decomposition through a sequence of abs-linear operations applied to intermediate quantities u. This can be done either directly on the pair of bounds ( u ˇ , u ^ ) or on their average u and their halved distance δ u = 1 2 ( u ˇ u ^ ) . In Section 3 we organize such sequences into an abs-linear form for f and then extend it to simultaneously yield the convex/concave decomposition. As a consequence of this analysis we get a strengthened version of the classical max min representation of piecewise linear functions, which reduces to the difference of two polyhedral parts in max- and min-form. In Section 4 we develop strict rules for propagating certain generalized gradient pairs ( g ˇ , g ^ ) of ( u ˇ , u ^ ) exploiting convexity and the cheap gradient principle [5]. In Section 5 we discuss the consequences for the DCA when using limiting gradients ( g ˇ , g ^ ) , solving the inner, linear optimization problem (LOP) exactly, and ensuring optimality via polyhedral reflection. In Section 6 we demonstrate the new results on the nonconvex and piecewise linear chained Rosenbrock version of Nesterov [6]. Section 7 contains a summary and preliminary conclusion with outlook. In the Appendix A we give the details of the necessary and sufficient optimality test from [7] in the present DC context.

2. Propagating Bounds and/or Radii

In Equation (3) we already assumed that doubling is done componentwise and that for a difference v = w u of DC functions w and u, one defines the convex and concave parts by
( w u ) ˇ = w ˇ u ^ and ( w u ) ^ = w ^ u ˇ ,
respectively. This yields in particular for the negation
( u ) ˇ = u ^ and ( u ) ^ = u ˇ .
For piecewise linear functions we need neither the square formula Equation (2) nor the more general decompositions for products. Therefore we will not insist on the sign conditions even though they would be also maintained automatically by Equation (4) as well as the natural linear rules for the convex and concave parts of the sum and the multiple of a DC function, namely
( w + u ) ˇ = ( w ˇ + u ˇ ) and ( w + u ) ^ = ( w ^ + u ^ ) , ( c u ) ˇ = c u ˇ and ( c u ) ^ = c u ^ if c 0 , ( c u ) ˇ = c u ^ and ( c u ) ^ = c u ˇ if c 0 .
However, the sign conditions would force one to decompose simple affine functions u ( x ) = a x + β  as
u ( x ) = max ( 0 , a x + β ) + min ( 0 , a x + β ) 1 2 ( u ˇ ( x ) + u ^ ( x ) ) ,
which does not seem such a good idea from a computational point of view.
The key observation for this paper is that as is well known (see e.g., [8]), one can propagate the absolute value operation according to the identity
| u | = max ( u , u ) = 1 2 max ( u ˇ + u ^ , u ˇ u ^ ) = max ( u ˇ , u ^ ) + 1 2 ( u ^ u ˇ ) | u | ˇ = 2 max ( u ˇ , u ^ ) and | u | ^ = u ^ u ˇ .
Here the equality in the second line can be verified by shifting the difference 1 2 ( u ^ u ˇ ) into the two arguments of the max. Again we see that when applying the absolute value operation to an already positive convex function u = 1 2 u ˇ 0 we get | u | ˇ = 2 u ˇ and | u | ^ = u ˇ so that the condition number grows from κ ( u ˇ , 0 ) = 1 to κ ( 2 u ˇ , u ˇ ) = 3 . In other words, we observe once more the danger that both component functions drift apart. This looks a bit like simultaneous growth of numerator and denominator in rational arithmetic, which can sometimes be limited through cancelations by common integer factors. It is currently not clear when and how a similar compactification of a given convex/concave decomposition can be achieved. The corresponding rule for the maximum is similarly easy derived, namely
max ( u , w ) = 1 2 max ( u ˇ + u ^ , w ˇ + w ^ ) = 1 2 max ( u ˇ w ^ , w ˇ u ^ ) + ( u ^ + w ^ ) .
When u and w as well as their decomposition are identical we arrive at the new decomposition u = max ( u , u ) = 1 2 ( ( u ˇ u ^ ) + 2 u ^ ) , which obviously represents again some deterioration in the conditioning.
While it was pointed out in [4] that the DC functions u = 1 2 ( u ˇ + u ^ ) themselves form an algebra, their decomposition pairs ( u ˇ , u ^ ) are not even an additive group, as only the zero ( 0 , 0 ) has a negative partner, i.e., an additive inverse. Naturally, the pairs ( u ˇ , u ^ ) form the Cartesian product between the convex cone of convex functions and its negative, i.e., the cone of concave functions. The DC functions are then the linear envelope of the two cones in some suitable space of locally Lipschitz continuous functions. It is not clear whether this interpretation helps in some way, and in any case we are here mainly concerned with piecewise linear functions.

Propagating the Center and Radius

Rather than propagating the pairs ( u ˇ , u ^ ) through an evaluation procedure as defined in [5] to calculate the function value f ( x ) at a given point x, it might be simpler and better for numerical stability to propagate the pair
u = 1 2 ( u ˇ + u ^ ) δ u = 1 2 ( u ˇ u ^ ) u ˇ = u + δ u u ^ = u δ u .
This representation resembles the so-called central form in interval arithmetic [9] and we will call therefore u the central value and δ u the radius. In other words, u is just the normal piecewise affine intermediate function and the δ u is a convex distance function to the hopefully close convex and concave part. Should the potential blow-up discussed above actually occur, this will only effect δ u but not the central value u itself. Moreover, at least theoretically one might decide to reduce δ u from time to time making sure of course that the corresponding u ˇ and u ^ as defined in Equation (7) stay convex and concave, respectively. The condition number now satisfies the bound
κ ( u + δ u , u δ u ) = sup x | u + δ u | + | u δ u | 2 | u | = sup x 1 2 | 1 + δ u u | + | 1 δ u u | 1 + sup x | δ u u | .
Recall here that all intermediate quantities u = u ( x ) are functions of the independent variable vector x R n . Naturally, we will normally only evaluate the intermediate pairs u and δ u at a few iterates of whatever numerical calculation one performs involving f so that we can only sample the ratio
ρ u ( x ) | δ u ( x ) / u ( x ) |
pointwise, where the denominator is hopefully nonzero. We will also refer to this ratio as the relative gap of the convex/concave decomposition at a certain evaluation point x. The arithmetic rules for propagating radii of the central form in central convex/concave arithmetic are quite simple.
Lemma 1 (Propagation rules for central form).
With c , d , x R two constants and an independent variable we have
v = c + d x δ v = 0 ρ v = 0 i f v 0 v = u ± w δ v = δ u + δ w ρ v | u | + | w | | u ± w | max ( ρ u , ρ w ) v = c u δ v = | c | δ u ρ v = ρ u i f c 0 v = | u | δ v = | u | + 2 δ u ρ v [ 1 , 1 + 2 ρ u ] .
Proof. 
The last rule follows from Equation (6) by
δ ( | u | ) = 1 2 | u | ˇ | u | ^ = max ( u ˇ , u ^ ) 1 2 ( u ^ u ˇ ) = max ( u ˇ δ u , u ^ δ u ) + 2 δ u = max ( u , u ) + 2 δ u = | u | + 2 δ u .
 ☐
The first equation in Equation (8) means that for all quantities u that are affine functions of the independent variables x the corresponding radius δ u is zero so that u ˇ = u = u ^ until we reach the first absolute value. Notice that δ v does indeed grow additively for the subtraction just like for the addition. By induction it follows from the rules above for an inner product that
δ j = 1 m c j u j = j = 1 m | c j | δ u j ,
where the c j R are assumed to be constants. As we can see from the bounds in Lemma 1 the relative gap can grow substantially whenever one performs an addition of values with opposite sign or applies the absolute value operation. In contrast to interval arithmetic on smooth functions one sees that the relative gap, though it may be zero or small initially immediately jumps above 1 when one hits the first absolute value operation. This is not really surprising since the best concave lower bound on u ( x ) = | x | itself is u ^ ( x ) = 0 so that δ u = | x | , u ˇ ( x ) = 2 | x | and thus ρ u ( x ) = 1 constantly. On the positive side one should notice that throughout we do not lose sight of the actual central values u ( x ) , which can be evaluated with full arithmetic precision. In any case we can think of neither ρ nor κ 1 + ρ as small numbers, but we must be content if they do not actually explode too rapidly. Therefore they will be monitored throughout our numerical experiments.
Again we see that the computational effort is almost exactly doubled. The radii can be treated as additional variables that occur only in linear operations and stay nonnegative throughout. Notice that in contrast to the (nonlinear) interval case we do not loose any accuracy by propagating the central form. It follows immediately by induction from Lemma 1 that any function evaluated by a evaluation procedure that comprises a finite sequence of
  • initializations to independent variables
  • multiplications by constants
  • additions or subtractions
  • absolute value applications
is piecewise affine and continuous. We will call these operations and the resulting evaluation procedure abs-linear. It is also easy to see that the absolute values | · | can be replaced by the maximum max ( · , · ) or the minimum min ( · , · ) or the positive part function max ( 0 , · ) or any combination of them, since they can all be mutually expressed in terms of each other and some affine operations. Conversely, it follows from the min-max representation established in [10] (Proposition  2.2.2) that any piecewise affine function f can be evaluated by such an evaluation procedure. Consequently, by applying the formulas Equations (4)–(6) one can propagate at the same time the convex and concave components for all intermediate quantities. Alternatively, one can propagate the centered form according to the rules given in Lemma 1. These rules are also piecewise affine so that we have a finite procedure for simultaneously evaluating u ˇ and u ^ or u and δ u as piecewise linear functions. The combined computation requires about 2–3 times as many arithmetic operations and twice as many memory accesses. Of course due to the interdependence of the two components it is not possible to evaluate just one of them without the other. As we will see the same is true for the generalized gradients to be discussed later in Section 4.

3. Forming and Extending the Abs-Linear Form

In practice all piecewise linear objectives can be evaluated by a sequence of abs-linear operations, possibly after min and max have been rewritten as
min ( u , w ) = 1 2 ( u + w | u w | ) and max ( u , w ) = 1 2 ( u + w + | u w | ) .
Our only restriction is that the number s of intermediate scalar quantities, say z i , is fixed, which is true for example in the max min representation. Then we can immediately cast the procedure in matrix-vector notation as follows:
Lemma 2 (Abs-Linear Form).
Any continuous piecewise affine function f : x R n y R can be represented by
z = c + Z x + M z + L | z | , y = d + a x + b z ,
where z R s , Z R s × n , M , L R s × s strictly lower triangular, d R , a R n , b R s and | z | denotes the componentwise modulus of the vector z.
It should be noted that the construction of this general abs-linear form requires no analysis or computation whatsoever. However, especially for our purpose of generating a reasonably tight DC decomposition, it is advantages to reduce the size of the abs-normal form by eliminating all intermediates z j with j < s for which | z j | never occurs on the right hand side. To this end we may simply substitute the expression of z j given in the j-th row in all places where z j itself occurs on the right hand side. The result is what we will call a reduced abs-normal form, where after renumbering, all remaining z j with j < s are switching variables in that | z j | occurs somewhere on the right hand side. In other words, all but the last column of the reduced, strictly lower triangular matrix L are nontrivial. Again, this reduction process is completely mechanical and does not require any nontrivial analysis, other than looking up which columns of the original L were zero. The resulting reduced system is smaller and probably denser, which might increase the computation effort for evaluating f itself. However, in view of Equation (9) we must expect that for the reduced form the radii will grow slower if we first accumulate linear coefficients and then take their absolute values. Hence we will assume in the remainder of this paper that the abs-normal form for our objective f of interest is reduced.
Based on the concept of abs-linearization introduced in [11], a slightly different version of a (reduced) abs-normal form was already proposed in [12]. Now in the present paper, both z and y depend directly on z via the matrix M and the vector b, but y does no longer depend directly on | z | . All forms can be easily transformed into each other by elementary modifications. The intermediate variables z i can be calculated successively for 1 i s by
z i = c i + Z i x + M i z + L i | z | ,
where Z i , M i and L i denote the ith rows of the corresponding matrix. By induction on i one sees immediately that they are piecewise affine functions z i = z i ( x ) , and we may define for each x the signature vector
σ ( x ) = ( sgn ( z i ( x ) ) ) i = 1 s { 1 , 0 , 1 } s .
Consequently we get the inverse images
P σ { x R n : sgn ( z ( x ) ) = σ } for σ { 1 , 0 , 1 } s ,
which are relatively open polyhedra that form collectively a disjoint decomposition of R n . The situation for the second example of Nesterov is depicted in Figure 3 in the penultimate section. There are six polyhedra of full dimension, seven polyhedra of co dimension 1 drawn in blue and two points, which are polyhedra of dimension 0. The point ( 0 , 1 ) with signature ( 0 , 1 , 0 ) is stationary and the point ( 1 , 1 ) with signature ( 1 , 0 , 0 ) is the minimizer as shown in [7]. The arrows indicate the path of our reflection version of the DCA method as described in Section 5.
When σ is definite, i.e., has no zero components, which we will denote by 0 σ , it follows from the continuity of z ( x ) that P σ has full dimension n unless it is empty. In degenerate situations this may also be true for indefinite σ but then the closure of P σ is equal to the extended closure
P σ ˜ { x R n : σ ( x ) σ ˜ } close ( P σ ˜ )
for some definite 0 σ ˜ σ . Here the (reflexive) partial ordering ≺ between the signature vectors satisfies the equivalence
σ ˚ σ σ ˚ i σ i σ i 2 for i = 1 s P σ ˚ P σ
as shown in [13]. One can easily check that for any σ σ ˚ there exists a unique signature
( σ σ ˚ ) i = σ i if σ ˚ i 0 σ i if σ ˚ i = 0 for i = 1 s
We call σ ˜ σ σ ˚ the reflection of σ at σ ˚ , which satisfies also σ ˜ σ ˚ and we have in fact P σ ˜ P σ = P σ ˚ . Hence the relation between σ and σ ˜ is symmetric in that also σ = σ ˜ σ ˚ . Therefore we will call ( σ , σ ˜ ) a complementary pair with respect to σ ˚ . In the very special case z i = x i for i = 1 n = s 1 the P σ are orthants and their reflections at the origin { 0 } = P 0 R n are their geometric opposites P σ ˜ with σ ˜ = σ . Here one can see immediately that all edges, i.e., one-dimensional polyhedra, have Cartesian signatures ± e i for i = 1 n and belong to P σ or P σ ˜ for any given σ . Notice that x ˚ is a local minimizer of a piecewise linear function if and and only if it is a local minimizer along all edges of nonsmoothness emanating form it. Consequently, optimality of f restricted to a complementary pair is equivalent to local optimality on R n , not only in this special case, but whenever the Linear Independence Kink Qualification (LIKQ) holds as introduced in [13] and defined in the Appendix A. This observation is the basis of the implicit optimality condition verified by our DCA variant Algorithm 1 through the use of reflections. The situation is depicted in Figure 3 where the signatures ( 1 , 1 , 1 ) and ( 1 , 1 , 1 ) as well as ( 1 , 1 , 1 ) and ( 1 , 1 , 1 ) form complementary pairs at ( 0 , 1 ) and ( 1 , 1 ) , respectively. At both reflection points there are four emanating edges, which all belong to one of the three polyhedra mentioned.
Applying the propagation rules from Lemma 1, one obtains with δ x = 0 R n the recursion
δ z 1 = δ ( c 1 + Z 1 x ) = 0 δ z i = ( | M i | + 2 | L i | ) δ z + | L i | | z | for i = 2 s ,
where the modulus is once more applied componentwise for vectors and matrices. Hence, we have again in matrix vector notation
δ z = ( | M | + 2 | L | ) δ z + | L | | z | ,
which yields for δ z the explicit expression
δ z = ( I | M | 2 | L | ) 1 | L | | z | = j = 0 ν ( | M | + 2 | L | ) j | L | | z | 0 .
Here, ν is the so-called switching depth of the abs-linear form of f, namely the largest ν N such that ( | M | + | L | ) ν 0 , which is always less than s due to the strict lower triangularity of M and L. The unit lower triangular ( I | M | 2 | L | ) is an M-matrix [14], and interestingly enough does not even depend on x but directly maps | z | = | z ( x ) | to δ z = δ z ( x ) . For the radius of the function itself, the propagation rules from Lemma 1 then yield
δ f ( x ) = δ y = | b | δ z 0 .
This nonnegativity implies the inclusion Equation (1) already mentioned in Section 1, i.e.:
Theorem 1 (Inclusion by convex/concave decomposition).
For any piecewise affine function f in abs-linear form, the construction defined in Section 2 yields a convex/concave inclusion
f ^ ( x ) f ( x ) 1 2 ( f ˇ ( x ) + f ^ ( x ) ) f ˇ ( x ) .
Moreover, the convex and the concave parts f ˇ ( x ) and f ^ ( x ) have exactly the same switching structure as f ( x ) in that they are affine on the same polyhedra P σ defined in (13).
Proof. 
Equations (16) and (17) ensure that δ f ( x ) is nonnegative at all x R n such that
f ^ ( x ) = f ( x ) δ f ( x ) f ( x ) f ( x ) + δ f ( x ) f ˇ ( x ) .
It follows from Equation (17) that the radii δ z i ( x ) are like the | z i ( x ) | piecewise linear with the only nonsmoothness arising through the switching variables z ( x ) themselves. Obviously this property is inherited by δ f ( x ) and the linear combinations f ˇ ( x ) = f ( x ) + δ f ( x ) and f ^ ( x ) = f ( x ) δ f ( x ) , which completes the proof. ☐
Combining Equations (16) and (18) with the abs-linear form of the piecewise affine function f and defining z ˜ = ( z , δ z ) R 2 s , one obtains for the calculation of f ˜ ( x ) y ˜ ( y , δ y ) the following abs-linear form
z ˜ = c ˜ + Z ˜ x + M ˜ z ˜ + L ˜ | z ˜ | ,
y ˜ = d ˜ + a ˜ x + b ˜ z ˜
with the vectors and matrices defined by
c ˜ = c 0 R 2 s , Z ˜ = Z 0 R 2 s × n , M ˜ = M 0 0 | M | + 2 | L | R 2 s × 2 s , L ˜ = L 0 | L | 0 R 2 s × 2 s , d ˜ = d 0 R 2 , a ˜ = a 0 R n × 2 , b ˜ = b 0 0 | b | R 2 s × 2 .
Then, Equations (19) and (20) yield
z δ z = c 0 + Z 0 x + M 0 0 | M | + 2 | L | z δ z + L 0 | L | 0 | z | | δ z | = c + Z x + M z + L | z | ( | M | + 2 | L | ) δ z + | L | | z | y δ y = d ˜ + a ˜ x + b ˜ z ˜ = d 0 + a x 0 + b 0 0 | b | z δ z = d + a x + b z | b | δ z ,
i.e., Equations (16) and (18). As can be seen, the matrices M ˜ and L ˜ have the required strictly lower triangular form. Furthermore, it is easy to check, that the switching depth of the abs-linear form of f carries over to the abs-linear form for f ˜ in that also ( | M ˜ | + | L ˜ | ) ν 0 = ( | M ˜ | + | L ˜ | ) ν + 1 . However, notice that this system is not reduced since the s radii are not switching variables, but globally nonnegative anyhow. We can now obtain explicit expressions for the central values, radii, and bounds for a given signature  σ .
Corollary 1 (Explicit representation of the centered form).
For any definite signature σ 0 and all x P σ we have with Σ = diag ( σ )
z σ ( x ) = ( I M L Σ ) 1 ( c + Z x ) a n d | z σ ( x ) | = Σ z σ ( x ) 0
δ z σ ( x ) = ( I | M | 2 | L | ) 1 | L | Σ ( I M L Σ ) 1 ( c + Z x ) 0
z σ = ( I M L Σ ) 1 Z σ f = a + b ( I M L Σ ) 1 Z
f σ = a + b + | b | ( I | M | 2 | L | ) 1 | L | Σ ( I M L Σ ) 1 Z
f ^ σ = a + b | b | ( I | M | 2 | L | ) 1 | L | Σ ( I M L Σ ) 1 Z ,
where the restrictions of the functions and their gradients to P σ are denoted by subscript σ. Notice that the gradients are constant on these open polyhedra.
Proof. 
Equations (21) and (23) follow directly from Equation (12), the abs-linear form (11) and the properties of Σ . Combining Equation (16) with  (21) yields Equation (22). Since f ˇ ( x ) = f ( x ) + δ f ( x ) and f ^ ( x ) = f ( x ) δ f ( x ) , Equations (24) and (25) follow from the representation in abs-linear form and Equation (23). ☐
As one can see the computation of the gradient f σ requires the solution of one unit upper triangular linear system and that of both f ˇ σ and f ^ σ one more. Naturally, upper triangular systems are solved by back substitution, which corresponds to the reverse mode of algorithmic differentiation as described in the following section. Hence, the complexity for calculating the gradients is exactly the same as that for calculating the functions, which can be obtained by one forward substitution for f σ and an extra one for δ f σ and thus f ˇ σ and f ^ σ . The given f σ , f ˇ σ and f ^ σ are proper gradients in the interior of the full dimensional domains P σ . For some or even many σ the inverse image P σ of the map x sgn ( z ( x ) ) may be empty, in which case the formulas in the corollary do not apply. Checking the nonemptiness of P σ for a given signature σ amounts to checking the consistency of a set of linear inequalities, which costs the same as solving an LOP and is thus nontrivial. Expressions for the generalized gradients at points in lower dimensional polyhedra are given in the following Section 4. There it is also not required that the abs-linear normal form has been reduced, but one may consider any given sequence of abs-linear operations.

The Two-Term Polyhedral Decomposition

It is well known ([15], Theorem 2.49) that all piecewise linear and globally convex or concave functions can be represented as the maximum or the minimum of a finite collection of affine functions, respectively. Hence, from the convex/concave decomposition we get the following drastic simplification of the classical min-max representation given, e.g., in [10].
Corollary 2 (Additive max/min decomposition of PL functions).
For every piecewise affine function f : R n R there exist k 0 affine functions α i + a i x for i = 1 k and l 0 affine functions β j + b j x for j = 1 l such that at all x R n
f ( x ) = max i = 1 k ( α i + a i x ) 1 2 f ˇ ( x ) + min j = 1 l ( β j + b j x ) 1 2 f ^ ( x )
where furthermore f ^ ( x ) f ( x ) f ˇ ( x ) .
The max-part of this representation is what is called a polyhedral function in the literature [15]. Since the min-part is correspondingly the negative of a polyhedral function we may also refer to Equation (26) as a DP decomposition, i.e., the difference of two polyhedral functions.
We are not aware of a publication that gives a practical procedure for computing such a collection of affine functions α i + a i x , i = 1 k , and β j + b j x , j = 1 l , for a given piecewise linear function f. Of course the critical question is in which form the function f is specified. Here as throughout our work we assume that it is given by a sequence of abs-linear operations. Then we can quite easily compute for each intermediate variable v representations of the form
v = i = 1 m ¯ max 1 j k i ( α i j + a i j x ) + i = 1 n ¯ min 1 j l i ( β i j + b i j x )
= max j i I i 1 i m ¯ i = 1 m ¯ ( α i j i + a i j i x ) + min j i J i 1 i n ¯ i = 1 n ¯ ( β i j i + b i j i x ) .
with index sets I i = { 1 , , k i } , 1 i m ¯ , and J i = { 1 , , l i } , 1 i n ¯ , since one has to consider all possibilities of selecting one affine function each from one of the m ¯ max and n ¯ min groups, respectively. Obviously, (28) involves i = 1 m k i and i = 1 n i affine function terms in contrast to the first representation (27) which contains just i = 1 m k i and i = 1 n i of them. Still the second version conforms to the classical representation of convex and concave piecewise linear functions, which yields the following result:
Corollary 3 (Explicit computation of the DP representation).
For any piecewise linear function given as abs-linear procedure one can explicitly compute the representation (26) by implementing the rules of Lemma 1.
Proof. 
We will consider the representations (27) from which (26) can be directly obtained in the form (28). Firstly, the independent variables x j are linear functions of themselves with gradient a = e j and inhomogeneity α = 0 . Then for multiplications by a constant c > 0 we have to scale all affine functions by c. Secondly, addition requires appending the expansions of the two summands to each other without any computation. Taking the negative requires switching the sign of all affine functions and interchanging the max and min group. Finally, to propagate through the absolute values we have to apply the rule (6), which means switching the signs in the min group, expressing it in terms of max and merging it with the existing max group. Here merging means pairwise joining each polyhedral term of the old max-group with each term in the switched min-group. Then the new min-group is the old one plus the old max-group with its sign switched. ☐
We see that taking the absolute value or, alternatively, maxima or minima generates the strongest growth in the number of polyhedral terms and their size. It seems clear that this representation is generally not very useful because the number of terms will likely blow up exponentially. This is not surprising because we will need one affine function for each element of the polyhedral decompositions of the domain of the max and min term. Typically, many of the affine terms will be redundant, i.e., could be removed without changing the values of the polyhedral terms. Unfortunately, identifying those already requires solving primal or dual linear programming problems, see, e.g., [16]. It seems highly doubtful that this would ever be worthwhile. Therefore, we will continue to advocate dealing with piecewise linear functions in a convenient procedural abs-linear representation.

4. Computation of Generalized Gradients and Constructive Oracle Paradigm

For optimization by variants of the DCA algorithm [17] one needs generalized gradients of the convex and the concave component. Normally, there are no strict rules for propagating generalized gradients through nonsmooth evaluation procedures. However, exactly this is simply assumed in the frequently invoked oracle paradigm, which states that at any point x R n the function value f ( x ) and an element g f ( x ) can be evaluated. We have argued in [18] that this is not at all a reasonable assumption.
On the other hand, it is well understood that for the convex operations: Positive scaling, addition, and taking the maximum the rules are strict and simple. Moreover, then the generalized gradient in the sense of Clarke f ˇ ( x ) R n is actually a subdifferential in that all its elements define supporting hyperplanes. Similarly f ^ ( x ) might be called a superdifferential in that the tangent planes bound the concave part from above.
In other words, we have at all x R n and for all increments Δ x
f ˇ ( x + Δ x ) f ˇ ( x ) + g ˇ Δ x if g ˇ f ˇ ( x )
and
f ^ ( x + Δ x ) f ^ ( x ) + g ^ Δ x if g ^ f ^ ( x ) ,
which imply for g ˇ f ˇ ( x ) and g ^ f ^ ( x ) that
f ^ ( x + Δ x ) + f ˇ ( x ) + g ˇ Δ x 2 f ( x + Δ x ) f ˇ ( x + Δ x ) + f ^ ( x ) + g ^ Δ x ,
where the lower bound on the left is a concave function and the upper bound is convex, both with respect to Δ x . Notice that the generalized superdifferential f ^ being the negative of the subdifferential of f ^ is also a convex set.
Now the key question is how we can calculate a suitable pair of generalized gradients ( g ˇ , g ^ ) f ˇ ( x ) × f ^ ( x ) . As we noted above the convex part and the negative of the concave part only undergo convex operations so that for v = c u
v ˇ = c u ˇ if c > 0 0 if c = 0 c u ^ if c < 0 and v ^ = c u ^ if c > 0 0 if c = 0 c u ˇ if c < 0
and for v = u + w
v ˇ = u ˇ + w ˇ and v ^ = u ^ + w ^ .
Finally, for v = | u | we find by Equation (6) that v ^ = u ^ u ˇ as well as
1 2 v ˇ = max ( u ˇ , u ^ ) = u ˇ if u > 0 conv u ˇ ( u ^ ) if u = 0 , u ^ if u < 0
where we have used that u = 1 2 ( u ˇ + u ^ ) in Equation (32). The sign of the arguments u of the absolute value function are of great importance, because they determine the switching structure. For this reason, we formulated the cases in terms of u rather than in the convex/concave components. The operator conv { · } denotes taking the convex hull or envelope of a given usually closed set. It is important to state that within an abs-linear representation the multipliers c will stay constant independent of the argument x, even if they were originally computed as partial derivatives by an abs-linearization process and thus subject to round-off error. In particular their sign will remain fixed throughout whatever algorithmic calculation we perform involving the piecewise linear function f. So, actually the case c = 0 could be eliminated by dropping this term completely and just initializing the left hand side v to zero.
Because we have set identities we can propagate generalized gradient pairs ( u ˇ , u ^ ) u ˇ × u ^ and perform the indicated algebraic operations on them, starting with the Cartesian basis vectors
x ˇ j = x ^ j = x j = e j since x ˇ j = x ^ j = x j for j = 1 n .
The result of this propagation is guaranteed to be an element of f ˇ × f ^ . Recall that in the merely Lipschitz continuous case generalized gradients cannot be propagated with certainty since for example the difference v = w u generates a proper inclusion v w u . In that vein we must emphasize that the average 1 2 ( f ˇ + f ^ ) need not be a generalized gradient of f = 1 2 ( f ˇ + f ^ ) as demonstrated by the possibility that f ^ = f ˇ algebraically but we happen to calculate different generalized gradients of f ˇ and f ^ at a particular point x. In fact, if one could show that f = 1 2 ( f ˇ + f ^ ) one would have verified the oracle paradigm, whose use we consider unjustified in practice. Instead, we can formulate another corollary for sufficiently piecewise smooth functions.
Definition 1.
For any d N , the set of functions f : R n R , y = f ( x ) , defined by an abs-normal form
z = F ( x , z , | z | ) , y = φ ( x , z ) ,
with F C d ( R n + s + s ) and φ C d ( R n + s ) , is denoted by C abs d ( R n ) .
Once more, this definition differs slightly from the one given in [7] in that y depends only on z and not on |z| in order to match the abs-linear form used here. Then one can show the following result:
Corollary 4 (Constructive Oracle Paradigm).
For any function f C abs 2 ( R n ) and a given point x there exist a convex polyhedral function Δ f ˇ ( x ; Δ x ) and a concave polyhedral function Δ f ^ ( x ; Δ x ) such that
f ( x + Δ x ) f ( x ) = 1 2 Δ f ˇ ( x ; Δ x ) + + Δ f ^ ( x ; Δ x ) + O ( Δ x 2 )
Moreover, both terms and their generalized gradients at Δ x = 0 or anywhere else can be computed with the same order of complexity as f itself.
Proof. 
In [11], we show that
f ( x + Δ x ) f ( x ) = Δ f ( x ; Δ x ) + O ( Δ x 2 ) ,
where Δ f ( x ; Δ x ) is a piecewise linearization of f developed at x and evaluated at Δ x . Applying the convex/concave decomposition of Theorem 1, one obtains immediately the assertion with a convex polyhedral function Δ f ˇ ( x ; Δ x ) and a concave polyhedral function Δ f ^ ( x ; Δ x ) evaluated at Δ x The complexity results follow from the propagation rules derived so far. ☐
We had hoped that it would be possible to use this approximate decomposition into polyhedral parts to construct at least locally an exact decomposition of a general function f C abs d ( R n ) . into a convex and compact part. The natural idea seems to add a sufficiently large quadratic term β ( Δ x 2 to
f ( x + Δ x ) f ( x ) 1 2 Δ f ^ ( x ; Δ x ) = 1 2 Δ f ˇ ( x ; Δ x ) + O ( Δ x 2 )
such that it would become convex. Then the same term could be subtracted from Δ f ^ ( x ; Δ x ) maintaining its concavity. Unfortunately, the following simple example shows that this is not possible.
Example 1 (Half pipe).
The function
f : R 2 R , f ( x 1 , x 2 ) = max ( x 2 2 max ( x 1 , 0 ) , 0 ) = x 2 2 i f x 1 0 x 2 2 x 1 i f 0 x 1 x 2 2 0 i f 0 x 2 2 x 1 ,
in the class C abs ( R n ) is certainly nonconvex as shown in Figure 1. As already observed in [19] this generally nonsmooth function is actually Fréchet differentiable at the origin x = 0 with a vanishing gradient f ( 0 ) = 0 Hence, we have f ( Δ x ) = O ( Δ x 2 and may simply choose constantly Δ f ˇ ( 0 ; Δ x ) 0 Δ f ^ ( 0 ; Δ x ) . However, neither by adding β Δ x 2 nor any other smooth function to f ( Δ x ) can we eliminate the downward facing kink along the vertical axis Δ x 1 = 0 . In fact, it is not clear whether this example has any DC decomposition at all.

Applying the Reverse Mode for Accumulating Generalized Gradients

Whenever gradients are propagated forward through a smooth evaluation procedure, i.e., for functions in C 2 ( R n ) , they are uniquely defined as affine combinations of each other, starting from Cartesian basis vectors for the components of x. Given only the coefficients of the affine combinations one can propagate corresponding adjoint values, or impact factors backwards, to obtain the gradient of a single dependent with respect to all independents at a small multiple of the operations needed to evaluate the dependent variable by itself. This cheap gradient result is a fundamental principle of computational mathematics, which is widely applied under various names, for example discrete adjoints, back propagation, and reverse mode differentiation. For a historical review see [20] and for a detailed description using similar notation to the current paper see our book [5]. For good reasons, there has been little attention to the reverse mode in the context of nonsmooth analysis, where one can at best obtain subgradients. The main obstacle is again that the forward propagation rules are only sharp when all elementary operations maintain convexity, which is by the way the only constructive way of verifying convexity for a given evaluation procedure. While general affine combinations and the absolute value are themselves convex functions, they do not maintain convexity when applied to a convex argument.
The last equation of Lemma 1 shows that one cannot directly propagate a subgradient of the convex radius functions δ u because there is a reference to v = | u | itself, which does not maintain convexity except when it is redundant due to its argument having a constant sign. However, it follows from the identity δ u = 1 2 ( u ˇ u ^ ) that for all intermediates u
u ˇ u ˇ u ^ u ^ 1 2 ( u ˇ u ^ ) δ u .
Hence one can get affine lower bounds of the radii, although one would probably prefer upper bounds to limit the discrepancy between the convex and concave parts. When v = | u | and u = 0 we may choose according to Equation (32) any convex combination
1 2 v ˇ = ( 1 μ ) u ˇ μ u ^ for 0 μ 1 .
It is tempting but not necessarily a good idea to always choose the weight μ equal to 1 2 for simplicity.
Before discussing the reasons for this at the end of this subsection, let us note that from the values of the constants c, the intermediate values u, and the chosen weights μ it is clear how the next generalized gradient pair ( v ˇ , v ^ ) is computed as a linear combination of the generalized gradients of the inputs for each operation, possibly with a switch in their roles. That means after only evaluating the function f itself, not even the bounds f ˇ and f ^ , we can compute a pair of generalized gradients in f ˇ × f ^ using the reverse mode of algorithmic differentiation, which goes back to at least [21] though not under that name. The complexity of this computation will be independent of the number of variables and relative to the complexity of the function f itself. All the operations are relatively benign, namely scaling by constants, interchanges and additions and subtractions. After all the reverse mode is just a reorganization of the linear algebra in the forward propagation of gradients. Hence, it appears that we can be comparatively optimistic regarding the numerical stability of this process.
To be specific we will indicate the (scalar) adjoint value of all intermediates u ˇ and u ^ as usual by u ˇ R and u ^ R . They are all initialized to zero except for either y ˇ = 1 or y ^ = 1 . Then at the end of the reverse sweep, the vectors ( x j ) j = 1 n represent either y ˇ or y ^ , respectively. For computational efficiency one may propagate both adjoint components simultaneously, so that one computes with sextuplets consisting of u ˇ , u ^ and their adjoints with respect to y ˇ and y ^ . In any case we have the following adjoint operations. For v = u + w
( w ˇ , w ^ ) += ( v ˇ , v ˇ ^ ) and ( u ˇ , u ^ ) += ( v ˇ , v ^ ) ,
for v = c u
( u ˇ , u ^ ) += c ( v ˇ , v ^ ) if c > 0 ( 0 , 0 ) if c = 0 c ( v ^ , v ˇ ) if c < 0 ,
and finally for v = | u |
( u ˇ , u ^ ) += ( 2 v ˇ v ^ , v ^ ) if u > 0 ( v ^ + 2 ( 1 μ ) v ˇ , v ^ 2 μ v ˇ ) if u = 0 ( v ^ , v ^ 2 v ˇ ) if u < 0 .
Of course, the update for the critical case u = 0 of the absolute value is just the convex combination for the two cases u > 0 and u < 0 weighted by μ. Due to round-off errors it is very unlikely that the critical case u = 0 ever occurs in floating point arithmetic. Once more, the sign of the arguments u of the absolute value function are of great importance, because they determine on which faces of the polyhedral functions f ˇ and f ^ the current argument x is located. In some situations one prefers a gradient that is limiting in that it actually occurs as a proper gradient on one of the adjacent smooth pieces. For example, if we had simply f ( x ) = v = | x | for x R and chose μ = 1 2 we would get v ˇ = 2 | x | , v ^ = 0 and find by Equation (34) that v ˇ = 2 ( 1 2 1 2 ) = 0 at x = x ˇ = x ^ = 0 . This is not a limiting gradient of v ˇ since v ˇ = [ 2 , 2 ] , whose interior contains the particular generalized gradient 0.

5. Exploiting the Convex/concave Decomposion for the DC Algorithm

In order to minimize the decomposed objective function f we may use the DCA algorithm [17] which is given in its basic form using our notation by
Choose x 0 R n For k = 0 , 1 , 2 ,         Calculate g k 1 2 f ^ ( x k )         Calculate x k + 1 1 2 f ˇ * ( g k )
where 1 2 f ˇ * denotes the Fenchel conjugate of 1 2 f ˇ . For a convex function h : R n R one has
w h * ( y ) w argmin x R n { h ( x ) y x } ,
see [15], Chapter 11. Hence, the classic DCA reduces in our Euclidean scenario to a simple recurrence
x k + 1 argmin x R n f ˇ ( x ) + g ^ k x for some g ^ k f ^ ( x k ) .
The objective function on the left hand side is a constantly shifted convex polyhedral upper bound on 2 f ( x ) since
f ˇ ( x ) + g ^ k x = 2 f ( x ) f ^ ( x ) g ^ k x 2 f ( x ) f ^ ( x k ) + g ^ k x k .
It follows from Equation (29) and x k + 1 being a minimizer that
f ( x k + 1 ) 1 2 f ˇ ( x k + 1 ) + f ^ ( x k ) + g ^ k ( x k + 1 x k ) 1 2 f ˇ ( x k ) + f ^ ( x k ) = f ( x k ) .
Now, since (36) is an LOP, an exact solution x k + 1 can be found in finitely many steps, for example by a variant of the Simplex method. Moreover, we can then assume that x k + 1 is one of finitely many vertex points of the epigraph of f ˇ . At these vertex points, f itself attains a finite number of bounded values. Provided f itself is bounded below, we can conclude that for any choice of the g ^ k f ^ σ ( k ) the resulting function values f ( x k ) can only be reduced finitely often so that f ( x k ) = f ( x k 1 ) and w.l.o.g. x k = x k 1 eventually. We then choose the next g ^ k = f ^ σ ( k ) with σ ( k ) = σ ( k 1 ) σ ( x k ) as the reflection of σ ( k 1 ) at σ ( x k ) as defined in (15). If then again f ( x k + 1 ) = f ( x k ) it follows from Corollary A2 that x k is a local minimizer of f and we may terminate the optimization run. Hence we obtain the DCA variant listed in Algorithm 1, which is guaranteed to reach local optimality under LIKQ. It is well defined even without this property and we conjecture that otherwise the final iterate is still a stationary point of f. The path of the algorithm on the example discussed in Section 5 is sketched in Figure 3. It reaches the stationary point ( 0 , 1 ) where σ = ( 0 , 1 , 0 ) from within the polyhedron with the signature ( 1 , 1 , 1 ) and then continues after the reflection ( 1 , 1 , 1 ) = ( 1 , 1 , 1 ) ( 0 , 1 , 0 ) . From within that polyhedron the inner loop reaches the point ( 1 , 1 ) with signature ( 1 , 0 , 0 ) , whose minimality is established after a search in the polyhedron P ( 1 , 1 , 1 ) .
If the function f ( x ) is unbounded below, so will be one of the inner convex problems and the convex minimizer should produce a ray of infinite descent instead of the next iterate x k + 1 . This exceptional scenario will not be explicitly considered in the remainder of the paper. The reflection operation is designed to facilitate further descent or establish local optimality. It is discussed in the context of general optimality conditions in the following subsection.
Algorithm1 Reflection DCA
Require: x 0 R n ,
1:
Set f 1 = and Evaluate f 0 = f ( x 0 )
2:
for k = 0 , 1 , do
3:
    if f k < f k 1 then                    ▹ Normal iteration with function reduction
4:
         Choose 0 σ σ ( x k )                  ▹ Here different heuristics may be applied
5:
         Compute g ^ k = f ^ σ                       ▹ Apply formula of Corollary 1
6:
    else                           ▹ The starting point was already optimal
7:
         Reflect σ ˜ = σ σ ( x k )                  ▹ The symbol ▹ is defined in Equation (15).
8:
         Update g ^ k = f ^ σ ˜
9:
    end if
10:
     Calculate x k + 1 argmin f ˇ ( x ) + g ^ k x x R n          ▹ Apply any LOP finite solver
11:
  Set f k + 1 = f ( x k + 1 )
12:
  if f k + 1 = f k = f k 1 then                      ▹ Local optimality established
13:
         Stop
14:
    end if
15:
endfor

5.1. Checking Optimality Conditions

Stationarity of x k happens when the convex function f ˇ ( x ) + g ^ k x is minimal at x k so that for all large k
0 f ˇ ( x k ) + g ^ k g ^ k f ^ ( x k ) ( f ˇ ( x k ) ) .
The nonemptiness condition on the right hand side is known as criticality of the DC decomposition at x k , which is necessary but not sufficient even for local optimality of f ( x ) at x k . To ensure the latter one has to verify that all g ^ k f ^ ( x k ) satisfy the criticality condition (38) so that
f ^ ( x k ) f ˇ ( x k ) L f ^ ( x k ) f ˇ ( x k ) .
The left inclusion is a well known local minimality condition [22], which is already sufficient in the piecewise linear case. The right inclusion is equivalent to the left one due to the convexity of f ˇ ( x k ) .
If f ˇ and f ^ were unrelated convex and concave polyhedral functions, one would normally consider it extremely unlikely that f ^ were nonsmooth at any one of the finitely many vertices of the polyhedral domain decomposition of f ˇ . For instance when f ^ is smooth at x k we find that f ^ ( x k ) = { g ^ k } is a singleton so that criticality according to Equation (38) is already sufficient for local minimality according to Equation (39). As we have seen in Theorem 1 the two parts have exactly the same switching structure. That means they are nonsmooth on the same skeleton of lower dimensional polyhedra. Hence, neither L f ˇ ( x k ) nor L f ^ ( x k ) will be singletons at minimizing vertices of the upper bound so that checking the validity of Equation (39) appears to be a combinatorial task at first sight.
However, provided the Linear Independence Kink Qualification (LIKQ) defined in [7] is satisfied at the candidate minimizer x k , the minimality can be tested with cubic complexity even in case of a dense abs-linear form. Moreover, if the test fails one can easily calculate a descent direction d. The details of the optimality test in our context including the calculation of a descent direction are given in the Appendix A. They differ slightly from the ones in [7]. Rather than applying the optimality test Proposition A1 explicitly, one can use its Corollary A2 stating that if x ˚ with σ ˚ = σ ( x ˚ ) is a local minimizer of the restriction of f to a polyhedron P σ with definite σ σ ˚ then it is a local minimizer of the unrestricted f if and only if it also minimizes the restriction of f to P σ ˜ with the reflection σ ˜ = σ σ ˚ . The latter condition must be true if x ˚ also minimizes f ( x ) + f ^ σ ˜ , which can be checked by solving that convex problem. If that test fails the optimization can continue.

5.2. Proximal Rather Than Global

By some authors the DCA algorithm has been credited with being able to reach global minimizers with a higher probability than other algorithms. There is really no justification for this optimism in the light of the following observation. Suppose the objective f ( x ) = 1 2 ( f ˇ ( x ) + f ^ ( x ) ) has an isolated local minimizer x * . Then there exists an ε > 0 such that the level set { x R n : f ( x ) f ( x * ) + ε } has a bounded connected component containing x * , say L ε . Now suppose DCA is started from any point x 0 L ε . Since f 0 ( x ) 1 2 ( f ˇ ( x ) + f ^ ( x 0 ) + g ^ ( x 0 ) ( x x 0 ) ) is by Equation (37) a convex upper bound on f ( x ) its level set { f 0 ( x ) f ( x 0 ) } will be contained in L ε . Hence any step from x 0 that reduces the upper bound f 0 ( x ) must stay in the same component, so there is absolutely no chance to move away from the catchment L ε of x 0 towards another local minimizer of f, whether global or not. In fact, by adding the convex term
1 2 f ^ ( x 0 ) + g ^ ( x 0 ) ( x x 0 ) f ^ ( x ) 0 ,
which vanishes at x 0 , to the actual objective f ( x ) one performs a kind of regularization, like in the proximal point method. This means the step is actually held back compared to a larger step that might be taken by a method that only requires the reduction of f ( x ) itself.
Hence we may interpret DCA as a proximal point method where the proximal term is defined as an affinely shifted negative of the concave part. Since in general the norm and the coefficient defining the proximal term may be quite hard to select, this way of defining it may make a lot of sense. However, it is certainly not global optimization. Notice that in this argument we have used neither the polyhedrality nor the inclusion property. So it applies to a general DC decomposition on Euclidean space. Another conclusion from the "holding back" observation is that it is probably not worthwhile to minimize the upper bound very carefully. One might rather readjust the shift g ^ x after a few or even just one iteration.

6. Nesterov’s Piecewise Linear Example

According to [6], Nesterov suggested three Rosenbrock-like test functions for nonsmooth optimization. One of them given by
f ( x ) = 1 4 | x 1 1 | + i = 1 n 1 x i + 1 2 | x i | + 1
is nonconvex and piecewise linear. It is shown in [6] that this function has 2 n 1 Clarke stationary points only one of which is a local and thus the global minimizer. Numerical studies showed that optimization algorithms tend to be trapped at one of the stationary points making it an interesting test problem. We have demonstrated in [23] that using an active signature strategy one can guarantee convergence to the unique minimizer from any starting point albeit using in the worst case 2 n iterations as all stationary points are visited. Let us first write the problem in the new abs-linear form.
Defining the s = 2 n switching variables
z i = F i ( x , | z | ) = x i for 1 i < n , z n = F n ( x , | z | ) = x 1 1 ,
and
z n + i = F n + i ( x , | z | ) = x i + 1 2 | z i | + 1 for 1 i < n , z s = 1 4 | z n | + i = 1 n 1 | z n + i |
the resulting objective function is then simply identical to y = f ( x ) = z s . With the vectors and matrices
c = ( 0 , 1 , e n 1 , 0 ) R ( n 1 ) + 1 + ( n 1 ) + 1 , Z = I n 1 0 I n 1 0 0 1 0 0 R s × ( n 1 ) + 1 , M = 0 , L = 0 0 0 0 0 0 0 0 2 I n 1 0 0 0 0 1 4 e n 1 0 R s × ( n 1 ) + 1 + ( n 1 ) + 1 , d = 0 R , a = 0 , b = ( 0 , , 0 , 1 ) R ( 2 n 1 ) + 1 ,
where Z and L have different row partitions, one obtains an abs-linear form (11) of f. Here, I k denotes the identity matrix of dimension k, e = ( 1 , , 1 ) R k the vector containing only ones and the symbol 0 pads with zeros to achieve the specified dimensions. One can easily check that | L | 2 0 = | L | 3 , hence this example has switching depth ν = 2 . The geometry of the situation is depicted in Figure 3, which was already briefly discussed in Section 3 and Section 5.
Since the corresponding extended abs-linear form for f ˜ = ( y , δ y ) does not provide any new insight we do not state it here. Directly in terms of the original equations we obtain for the radii
δ z i = 0 for 1 i n , δ z n + i = 2 | z i | = 2 | x i | for 1 i < n
and
δ f = δ z s = 1 4 | z n | + i = 1 n 1 ( | z n + i | + 2 δ z n + i ) = 1 4 | x 1 1 | + i = 1 n 1 ( | x i + 1 2 | x i | + 1 | + 4 | x i | ) .
Thus, from Equation (7) we get the convex and concave part explicitly as
z ˇ i = z i = z ^ i for 1 i n , z ˇ n + i = x i + 1 + 1 z ^ n + i = x i + 1 4 | z i | + 1 = x i + 1 4 | x i | + 1 for 1 i < n
and most importantly
f ˇ = z s + δ z s = 1 2 | x 1 1 | + 2 i = 1 n 1 x i + 1 2 | x i | + 1 + 2 | x i | f ^ = z s δ z s = 4 i = 1 n 1 | x i | .
Clearly f ^ is a concave function and to check the convexity of f ˇ we note that
x i + 1 2 | x i | + 1 + 2 | x i | = 2 | x i | 1 x i + 1 + 2 | x i | 1 x i + 1 + x i + 1 + 1 = 1 + x i + 1 + 2 max 0 , 2 | x i | x i + 1 1 .
The last expression is the sum of an affine function and the positive part of the sum of the absolute value and an affine function, which must therefore also be convex. The corresponding term in Equation (42) is the same with the convex function 2 | x i | added, so that δ f is also convex in agreement with the general theory. Finally, one verifies easily that
f ^ f = 1 2 ( f ˇ + f ^ ) f ˇ ,
which is the whole idea of the decomposition. It would seem that the automatic decomposition by propagation through the abs-linear procedure yields a rather tight result. The function f as well as the lower and upper bound given by the convex/concave decomposition are illustrated on the left hand side of Figure 2. Notice that the switching structure is indeed identical for all three as stated in Theorem 1. On the right hand side of Figure 2, the difference 2 δ f between the upper, convex and lower, concave bound is shown, which is indeed convex.
It is worthwhile to look at the condition number of the decomposition, namely we get the following trivial bound
κ ( f ˇ , f ^ ) = sup x R n 1 2 | x 1 1 | + 2 i = 1 n 1 x i + 1 2 | x i | + 1 + 4 | x i ) 1 2 | x 1 1 | + 2 i = 1 n 1 | x i + 1 2 | x i | + 1 | = 1 + sup x R n 8 i = 1 n 1 | x i | 1 4 | x 1 1 | + 2 i = 1 n 1 | x i + 1 2 | x i | + 1 | = .
The disappointing right hand side value follows from the fact that at the well known unique global optimizer x * = ( 1 , 1 , , 1 ) R n the numerator is zero and the denominator positive. However, elsewhere, we can bound the conditioning as follows.
Lemma 3.
In case of the example (40) there is a constant c R such that
κ ( f ˇ ( x ) , f ^ ( x ) ) 1 + c min ( x x * , 3 ) .
Proof. 
Since the denominator is piecewise linear and vanishes only at the minimizer x * there must be a constant c 0 > 0 such that for x x * 3
8 i = 1 n 1 | x i | 1 4 | x 1 1 | + 2 i = 1 n 1 | x i + 1 2 | x i | + 1 | 8 i = 1 n 1 | x i | c 0 x x * 8 ( n 1 ) x c 0 x x * 32 ( n 1 ) c 0 x x * ,
which takes the value 32 ( n 1 ) / ( 3 c 0 ) on the boundary. On the other hand we get for x 2 and thus in particular x x * 3
8 i = 1 n 1 | x i | 1 4 | x 1 1 | + 2 i = 1 n 1 | x i + 1 2 | x i | + 1 | 4 ( n 1 ) x max 1 i < n | 2 | x i | x i + 1 1 | 4 ( n 1 ) 2 1 1 / 2 8 ( n 1 ) .
Assuming without loss of generality that c 0 4 / 3 we can combine the two bounds to obtain the assertion with c 32 ( n 1 ) / c 0 . ☐
Hence, we see the condition number κ ( f ˇ ( x ) , f ^ ( x ) ) is nicely bounded and the decomposition should work as long as our optimization algorithm has not yet reached its goal x * . It is verified in the companion article [24], that the DCA exploiting the observations made in this paper reaches the global minimizer in finitely many steps. It was already shown in [7] that the LIKQ condition is satisfied everywhere and that the optimality test singles out the unique minimizer correctly. In Figure 3, the arrows indicate the path of our reflection version of the DCA method as described in Section 5.

7. Summary, Conclusions and Outlook

In this paper the following new results were achieved
  • For every piecewise linear function f given as an abs-linear evaluation procedure, rules for simultaneously evaluating its representation as the average of a concave lower bound f ^ and a convex upper bound f ˇ are derived.
  • The two bounds can be constructively expressed as a single maximum and minimum of affine functions, which drastically simplifies the classical min max representation. Due to its likely combinatorial complexity we do not recommend this form for practical calculations.
  • For the two bounds f ˇ and f ^ , generalized gradients g ˇ and g ^ can be propagated forward or reverse through the convex or concave operations that define them. The gradients are not unique but guaranteed to yield supporting hyperplanes and thus provide a verified version of the oracle paradigm.
  • The DCA algorithm can be implemented such that a local minimizer is reached in finitely many iterations, provided the Linear Independence Kink Qualification (LIKQ) is satisfied. It is conjectured that without this assumption the algorithm still converges in finitely many steps to a Clarke stationary point. Details on this can be found in the companion paper [24].
These results are illustrated on the piecewise linear Rosenbrock variant of Nesterov.
On a theoretical level it would be gratifying and possibly provide additional insight, to prove the result of Corollary A3 directly using the explicit representations of the generalized differentials of the convex and concave part given in Corollary 1. Moreover, it remains to be explored what happens when LIKQ is not satisfied. We have conjectured in [25] that just verifying the weaker Mangasarian Fromovitz Kink Qualification (MFKQ) represents an NP hard task. Possibly, there are other weaker conditions that can be cheaply verified and facilitate the testing for at least local optimality.
Global optimality can be characterized theoretically in terms of ε subgradients, albeit with ε arbitrarily large [26]. There is the possibility that the alternative definition of ε-gradients given in [18] might allow one to constructively check for global optimality. It does not really seem clear how these global optimality conditions can be used to derive corresponding algorithms.
The implementation of the DCA algorithm can be optimized in various ways. Notice that for applying the Simplex method in standard form, one could use for the representation as DC function the max-part in the more economical representation Equation (27) introducing m ¯ additional variables, rather than the potentially combinatorial Equation (28) to assemble the constraint matrix. In any case it seems doubtful that solving each sub problem to completion is a good idea, especially as the resulting step in the outer iteration is probably much too small anyhow. Therefore, the generalized gradient of the concave part, which defines the inner problem, should probably be updated much more frequently. Moreover, the inner solver might be an SQOP type active signature method or a matrix free gradient method with momentum term, as is used in machine learning, notwithstanding the nonsmoothness of the objective. Various options in that range will be discussed and tested in the companion article [24].
Finally, one should always keep in mind that the task of minimizing a piecewise linear function will most likely occur as an inner problem in the optimization of a piecewise smooth and nonlinear function. As we have shown in [27] the local piecewise linear model problem can be obtained easily by a slight generalization of automatic or algorithmic differentiation, e.g., ADOL-C [28] and Tapenade [29].

Author Contributions

Conceptualization, A.G. and A.W.; methodology, A.G. and A.W.; writing—original draft preparation; writing—review and editing, A.G. and A.W. All authors have read and agreed to the published version of the manuscript.

Funding

We acknowledge support by the German Research Foundation (DFG) and the Open Access Publication Fund of Humboldt-Universität zu Berlin.

Acknowledgments

We thank Napsu Karmitsa and Sona Taheri for inviting us to participate in this special issue in honor of Adil M. Bagirov. We also thank the three anonymous referees, who asked for various corrections and clarifications, which made the paper much more self-contained and readable.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Polynomial Optimality Test Based on Abs-Linear Form

As illustrated for the Nesterov test function, it may be advantageous to use intermediate variables z i that are not arguments of the absolute value themselves. For simplicity, we assume that these switching variables that do not impose nonsmoothness are located in the last components of z and that only the s ˜ s components z 1 , z s ˜ are arguments of the absolute value. Let us abbreviate the current iterate x k with x ˚ x k and denote the corresponding switching vector by z ˚ = z ( x ˚ ) , the signature vector σ ˚ = sgn ( z ˚ ) and the active index set by α { i s ˜ : σ ˚ i = 0 } with cardinality m | α | s ˜ . Consequently, there are exactly 2 m definite signatures by σ σ ˚ and the same number of limiting gradients for the three generalized differentials f ˇ , f ^ , and f .
For all x P σ ˚ , the signature σ ˚ is constant and we can use Corollary 1 to define the smooth function
z σ ˚ ( x ) = ( I M L Σ ˚ ) 1 ( c + Z x ) = c ˚ + Z ˚ x ,
where we have pulled out the unit lower triangular factor ( I M L Σ ˚ ) such that
Z ˚ = ( I M L Σ ˚ ) 1 Z and c ˚ = ( I M L Σ ˚ ) 1 c .
For x x ˚ to be contained in the extended closure P σ ˚ as defined in Equation (14), it must satisfy the m linear equations
P α z ( x ) = 0 R m for P α = ( e i ) i α R m × s ˜
with e i denoting the ith unit vector in R s ˜ . Thus it is necessary and sufficient for P σ ˚ to be a polyhedron of dimension n m that the Jacobian P α Z ˚ R m × n has full row rank m. This rank condition was introduced as LIKQ in [7] and obviously requires that no more than n switches are active at x ˚ . As discussed in [7], for the point x ˚ to be a local minimizer of f it is necessary that it solves the trunk problem
min a x + b z s . t . | Σ ˚ | z c ˚ Z ˚ x = 0 .
Here | Σ ˚ | R s ˜ × s ˜ is the projection onto the s ˜ m vector components whose indices do not belong to α so the equality constraint combines (A1) and the constraint P α z = 0 . Now we get from KKT theory or equivalently LOP duality that x ˚ is a minimizer on P α if and only if for some Lagrange multiplier vector λ R s ˜
a = λ Z ˚ and b = λ | Σ ˚ | .
Since I = | Σ ˚ | + P α P α we derive that
λ ( I | Σ ˚ | ) Z ˚ = λ α P α Z ˚ = a b Z ˚ .
where λ α P α λ . This is a generally overdetermined system of n equations in the m components of λ α . If it is solvable the full multiplier vector λ = P α λ α + | Σ ˚ | b is immediately available. Because of the assumed full rank of the Jacobian P α Z ˚ we have m n , and if x ˚ is a vertex in that m = n the tangential stationarity condition (A3) is automatically satisfied.
Now it is necessary and sufficient for local minimality that x ˚ is also a minimizer of f on all polyhedra P σ with definite σ σ ˚ . Any such σ σ ˚ can be written as σ = σ ˚ + γ with γ { 1 , 0 , 1 } s ˜ structurally orthogonal to σ ˚ such that for Γ = diag ( γ ) we have the matrix equations
Σ = Σ ˚ + Γ and Σ ˚ Γ = 0 = | Σ ˚ | Γ .
Then we can express the z ( x ) = z σ ( x ) for x P σ as
z σ ( x ) = z σ ˚ + γ ( x ) = ( I M L Σ ˚ L Γ ) 1 ( c + Z x ) = ( I L ˚ Γ ) 1 ( c ˚ + Z ˚ x ) ,
with L ˚ ( I M L Σ ˚ ) 1 L . Now x ˚ must be the minimizer of f on P σ , i.e., solve the problem
min a x + b z s . t . ( I L ˚ Γ ) z = c ˚ + Z ˚ x , P α Γ z 0 R m .
Notice that the inequalities are only imposed on the sign constraints that are active at x ˚ since the strict inequalities are maintained in a neighborhood of x ˚ due to the continuity of z ( x ) . Then we get again from KKT theory or equivalently LOP duality that still a = λ Z ˚ and for a second multiplier vector 0 μ R m the equalities
a = λ Z ˚ and b = λ ( I L ˚ Γ ) + μ P α Γ .
Multiplying from the right by the projection | Σ ˚ | we find that the conditions (A2) and (A3) must still hold so that λ remains exactly the same. Moreover, multiplying from the right by Γ P α we get with P α P α = I m and Γ Γ = P α P α after some rearrangement the inequality
( λ b ) Γ P α = λ L ˚ P α μ λ L ˚ P α .
Now the key observation is that this condition is linear in Γ and is strongest for the choice γ i = sgn ( λ i b i ) for i α yielding the inequalities
| λ i b i | e i L ˚ λ for i α .
In other words, x ˚ is a solution of the branch problems (A4) if and only if it is for the worst case where γ i = sgn ( λ i b i ) for i α . When coincidentally λ i = b i we can define γ i arbitrarily. Note that the complementarity condition μ P α z ( x ˚ ) = 0 associated with Equation (A4) is automatically satisfied at x ˚ for any μ, since P α z ˚ = 0 by definition of the active index set α. These observations yield immediately:
Proposition A1 (Necessary and sufficient minimality condition).
Assume LIKQ holds in that P α Z ˚ has full row rank m = | α | . Then the point x ˚ is a local minimizer of f if and only if we have tangential stationarity in that a + Z ˚ b belongs to the range of Z ˚ P α and normal growth holds in that | P α ( λ b ) | P α L ˚ λ .
The verification that LIKQ holds and subsequently the test whether tangential stationarity is satisfied can be based on a QR decomposition of the active Jacobian P α Z ˚ R m × n . The main expense here is the calculation of Z ˚ itself, which requires one forward substitution on ( I M L Σ ˚ ) for each of n columns of Z and hence at most n s 2 / 2 fused multiply adds. Very likely this effort will already be made by any kind of active set method for reaching the candidate point x ˚ . Once the multiplier vector λ is obtained the remaining test (A7) for normal growth is almost for free so that we have a polynomial minimality criterion provided LIKQ holds. Otherwise one may assume a weaker generalization of the Mangasarian Fromovitz constrained qualification called MFKQ in [25]. However, we have conjectured in [19] that verifying MFKQ is probably already NP-hard.
Corollary A1 (Descent direction in the nonoptimal case).
Suppose that LIKQ holds. If tangential stationarity is violated there exits some direction d R n such that P α Z ˚ d = 0 but ( a + b Z ˚ ) d < 0 , which implies descent in that f ( x ˚ + τ d ) < f ( x ˚ ) for τ 0 . If tangential stationarity holds but normal growth fails there exists at least one i α with | λ i b i | > e i L ˚ λ . Defining γ = sgn ( λ i b i ) e i R s ˜ , any d satisfying P α ( I L ˚ Γ ) 1 Z ˚ d = P α γ is a descent direction.
Proof. 
In the first case it is clear that x ˚ + τ d P σ ˚ for τ 0 since the components of z ( x ˚ + τ d ) with indices in α stay zero and the others vary only slightly. Then the directional derivative of f ( . ) at x ˚ in direction τ d is given by
τ a d + τ b Z ˚ d = τ ( a d + b Z ˚ d ) < 0 ,
which proves the first assertion. Otherwise, λ is well defined and we can choose i α with | λ i b i | > e i L ˚ λ . Setting γ = γ i e i with γ i = sgn ( λ i b i ) e i , one obtains for d with P α ( I L ˚ Γ ) 1 Z ˚ d = γ that x ˚ + τ d P σ ˚ + γ for τ 0 . On that polyhedron the Lagrange multiplier vector μ is also well defined by Equation (A6) but we have
μ i = e i L ˚ λ ( λ i b i ) γ i = e i L ˚ λ | λ i b i | < 0 .
Then we get the directional derivative of f ( . ) at x ˚ in direction τ d
τ a d + τ b ( I L ˚ Γ ) 1 Z ˚ d = τ ( λ Z ˚ d + λ Z ˚ d + μ P α Γ ( I L ˚ Γ ) 1 Z ˚ d ) = τ μ i γ i 2 < 0 ,
where we have used identity (A5). Hence we have again descent, which completes the proof.  ☐
Corollary A2 (Optimality via Reflection).
Suppose an x ˚ where LIKQ holds has been reached by minimizing f ( x ) + g ^ x with g ^ = f ^ σ for 0 σ σ ˚ . Then x ˚ is a local minimizer of f on R n if and only if it is also a minimizer of f ˇ ( x ) + f ^ σ ˜ x with σ ˜ = σ σ ˚ as defined in (15).
Proof. 
By assumption x ˚ solves one of the branch problems of f itself. Hence we must have tangential stationarity (A5) with the corresponding Γ = diag ( γ ) for γ = σ σ ˚ . Since σ ˜ σ ˚ = γ we conclude from (A6) that
( λ b ) Γ P α λ L ˚ P α ( λ b ) ( Γ ) P α = ( λ b ) Γ P α
which implies that
( λ b ) P α = ( λ b ) Γ P α λ L ˚ P α .
Hence both tangential stationarity and normal growth are satisfied, which completes the proof by Proposition A1 as the converse implication is trivial. ☐
The key conclusion is that if an x ˚ is the solution of two complementary convex problems it must be locally optimal in the full dimensional space R n . Hence one can establish local optimality just using the preferred convex solver. If this test fails one naturally obtains descent to function values below f ( x ˚ ) until eventually a local minimizer is found.

Equivalence to DC Optimality Condition

Using the explicit expressions given in Lemma 1 we find that (see [18])
L f ( x ˚ ) = 0 = γ σ ˚ a + b ( I L ˚ Γ ) 1 Z ˚ ,
where γ ranges over all complements of σ ˚ such that σ ˚ + γ { 1 , 1 } s is definite. Similarly we obtain with
b ˜ | b | ( I | M | 2 | L | ) 1 | L | 0 R s
the limiting differentials of the convex and the concave part as
L f ˇ ( x ˚ ) = 0 = γ σ ˚ a + ( b + b ˜ Σ ˚ + b ˜ Γ ) ( I L ˚ Γ ) 1 Z ˚ ,
L f ^ ( x ˚ ) = 0 = γ σ ˚ a + ( b b ˜ Σ ˚ b ˜ Γ ) ( I L ˚ Γ ) 1 Z ˚ .
Hence we have an explicit representation for the limiting gradients of f as well as its convex and concave part f ˇ and f ^ at x ˚ . It is easy to see that the minimality condition (A5) requires a to be in the range of Z ˚ so that we have again a = λ Z ˚ yielding
L f ˇ ( x ˚ ) = 0 = γ σ ˚ ( b λ + λ L ˚ Γ + b ˜ Σ ˚ + b ˜ Γ ) ( I L ˚ Γ ) 1 Z ˚ ,
L f ^ ( x ˚ ) = 0 = γ σ ˚ ( b λ + λ L ˚ Γ b ˜ Σ ˚ b ˜ Γ ) ( I L ˚ Γ ) 1 Z ˚ .
We had hoped to be able to derive directly from these expressions that normal growth implies the condition (39), but we have so far not been able to do so. However, we can indirectly derive the following equivalence.
Corollary A3 (First order minimality condition).
Under LIKQ the limiting differential L f ^ ( x ˚ ) is contained in the convex hull of L f ˇ ( x ˚ ) if and only if tangential stationarity and normal growth condition hold according to Proposition A1.

References

  1. Joki, K.; Bagirov, A.; Karmitsa, N.; Mäkelä, M. A proximal bundle method for nonsmooth DC optimization utilizing nonconvex cutting planes. J. Glob. Optim. 2017, 68, 501–535. [Google Scholar] [CrossRef]
  2. Tuy, H. DC optimization: Theory, methods and algorithms. In Handbook of Global Optimization; Springer: Boston, MA, USA, 1995; pp. 149–216. [Google Scholar]
  3. Rump, S. Fast and parallel interval arithmetic. BIT 1999, 39, 534–554. [Google Scholar] [CrossRef]
  4. Bačák, M.; Borwein, J. On difference convexity of locally Lipschitz functions. Optimization 2011, 60, 961–978. [Google Scholar]
  5. Griewank, A.; Walther, A. Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation; Society for Industrial and Applied Mathematics: Philadelphia, PA, USA, 2008. [Google Scholar]
  6. Gürbüzbalaban, M.; Overton, M. On Nesterov’s nonsmooth Chebyshev-Rosenbrock functions. Nonlinear Anal. Theory Methods Appl. 2012, 75, 1282–1289. [Google Scholar]
  7. Griewank, A.; Walther, A. First and second order optimality conditions for piecewise smooth objective functions. Optim. Methods Softw. 2016, 31, 904–930. [Google Scholar] [CrossRef]
  8. Strekalovsky, A. Local Search for Nonsmooth DC Optimization with DC Equality and Inequality Constraints. In Numerical Nonsmooth Optimization. State of the Art Algorithms; Springer Nature Switzerland AG: Cham, Switzerland, 2020; pp. 229–261. [Google Scholar]
  9. Hansen, E. (Ed.) The centred form. In Topics in Interval Analysis; Oxford University Press: Oxford, UK, 1969; pp. 102–105. [Google Scholar]
  10. Scholtes, S. Introduction to Piecewise Differentiable Functions; Springer: New York, NY, USA, 2012. [Google Scholar]
  11. Griewank, A. On Stable Piecewise Linearization and Generalized Algorithmic Differentiation. Optim. Methods Softw. 2013, 28, 1139–1178. [Google Scholar] [CrossRef]
  12. Griewank, A.; Bernt, J.U.; Radons, M.; Streubel, T. Solving piecewise linear equations in abs-normal form. Linear Algebra Appl. 2015, 471, 500–530. [Google Scholar] [CrossRef]
  13. Griewank, A.; Walther, A.; Fiege, S.; Bosse, T. On Lipschitz optimization based on gray-box piecewise linearization. Math. Program. Ser. A 2016, 158, 383–415. [Google Scholar] [CrossRef]
  14. Golub, G.; Van Loan, C. Matrix Computations, 4th ed.; Johns Hopkins University Press: Baltimore, MD, USA, 2013. [Google Scholar]
  15. Rockafellar, R.; Wets, R.B. Variational Analysis; Springer: Berlin/Heidelberg, Germany, 1998. [Google Scholar]
  16. Fukuda, K.; Gärtner, B.; Szedlák, M. Combinatorial redundancy detection. Ann. Oper. Res. 2018, 265, 47–65. [Google Scholar] [CrossRef] [Green Version]
  17. Le Thi, H.; Pham Dinh, T. DC programming and DCA: Thirty years of developments. Math. Program. Ser. B 2018, 169, 5–68. [Google Scholar] [CrossRef]
  18. Griewank, A.; Walther, A. Beyond the Oracle: Opportunities of Piecewise Differentiation. In Numerical Nonsmooth Optimization. State of the Art Algorithms; Springer: Cham, Switzerland, 2020; pp. 331–361. [Google Scholar]
  19. Walther, A.; Griewank, A. Characterizing and testing subdifferential regularity in piecewise smooth optimization. SIAM J. Optim. 2019, 29, 1473–1501. [Google Scholar] [CrossRef]
  20. Griewank, A. Who invented the reverse mode of differentiation? Doc. Math. 2012, 389–400. [Google Scholar]
  21. Linnainmaa, S. Taylor expansion of the accumulated rounding error. BIT 1976, 16, 146–160. [Google Scholar] [CrossRef]
  22. Sun, W.; Sampaio, R.; Candido, M. Proximal point algorithm for minimization of DC function. J. Comput. Math. 2003, 21, 451–462. [Google Scholar]
  23. Griewank, A.; Walther, A. Finite convergence of an active signature method to local minima of piecewise linear functions. Optim. Methods Softw. 2019, 34, 1035–1055. [Google Scholar] [CrossRef]
  24. Griewank, A.; Walther, A. The True Steepest Descent Methods Revisited; Technical Report; Humboldt-Universität zu Berlin: Berlin, Germany, 2020. [Google Scholar]
  25. Griewank, A.; Walther, A. Relaxing kink qualifications and proving convergence rates in piecewise smooth optimization. SIAM J. Optim. 2019, 29, 262–289. [Google Scholar] [CrossRef] [Green Version]
  26. Niu, Y. Programmation DC & DCA en Optimisation Combinatoire et Optimisation Polynomiale via les Techniques de SDP. Ph.D. Thesis, INSA Rouen, Rouen, France, 2010. [Google Scholar]
  27. Fiege, S.; Walther, A.; Kulshreshtha, K.; Griewank, A. Algorithmic differentiation for piecewise smooth functions: A case study for robust optimization. Optim. Methods Softw. 2018, 33, 1073–1088. [Google Scholar] [CrossRef]
  28. Walther, A.; Griewank, A. Getting Started with ADOL-C. In Combinatorial Scientific Computing; Chapman & Hall/CRC Computational Science Series; CRC Press: Boca Raton, FL, USA, 2012; pp. 181–202. [Google Scholar]
  29. Hascoët, L.; Pascual, V. The Tapenade Automatic Differentiation tool: Principles, Model, and Specification. ACM Trans. Math. Softw. 2013, 39, 20:1–20:43. [Google Scholar]
Figure 1. Half pipe example as defined in Equation (33).
Figure 1. Half pipe example as defined in Equation (33).
Algorithms 13 00166 g001
Figure 2. Nesterov–Rosenbrock test function polyhedral inclusion for n = 2 .
Figure 2. Nesterov–Rosenbrock test function polyhedral inclusion for n = 2 .
Algorithms 13 00166 g002
Figure 3. Signatures and reflection-based DCA for Nesterov–Rosenbrock variant (40) with n = 2 .
Figure 3. Signatures and reflection-based DCA for Nesterov–Rosenbrock variant (40) with n = 2 .
Algorithms 13 00166 g003

Share and Cite

MDPI and ACS Style

Griewank, A.; Walther, A. Polyhedral DC Decomposition and DCA Optimization of Piecewise Linear Functions. Algorithms 2020, 13, 166. https://0-doi-org.brum.beds.ac.uk/10.3390/a13070166

AMA Style

Griewank A, Walther A. Polyhedral DC Decomposition and DCA Optimization of Piecewise Linear Functions. Algorithms. 2020; 13(7):166. https://0-doi-org.brum.beds.ac.uk/10.3390/a13070166

Chicago/Turabian Style

Griewank, Andreas, and Andrea Walther. 2020. "Polyhedral DC Decomposition and DCA Optimization of Piecewise Linear Functions" Algorithms 13, no. 7: 166. https://0-doi-org.brum.beds.ac.uk/10.3390/a13070166

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop