Between Piecewise Smoothness and Linear Complementarity Andreas - - PowerPoint PPT Presentation

Between Piecewise Smoothness and Linear Complementarity

Andreas Griewank, with thanks to

• A. Walther, T. Bosse, N. Strogies, S. Fiege, F. Kerkoff, J.U. Bernt,
• M. Radons,T. Streubel, R. Hasenfelder, P. Boeck, B. Lenser, . . .

Institute for Applied Mathematics, Humboldt Universit¨ at zu Berlin, Germany griewank@math.hu-berlin.de 6th International Conference on Complementarity Problems August 4-8, 2014, Berlin.

August 8, 2014

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 1 / 43

1 Observations/Opinions of a Johnny Come Lately 2 Piecewise Linearization/Differentiation 3 Representation of PL functions in abs-normal form 4 Computation of conical Jacobians and gradients 5 Solving PL systems of Equation and LCPs 6 (Un)constrained optimization by successive PL 7 Integration of Lipschitzian dynamics

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 2 / 43

Observations/Opinions of a Johnny Come Lately

Why not smoothen/regularize/parametrize?

• If scale too small it makes no algorithmic difference.

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 3 / 43

Observations/Opinions of a Johnny Come Lately

Why not smoothen/regularize/parametrize?

• If scale too small it makes no algorithmic difference.
• If scale too large the problem is significantly changed.

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 3 / 43

Observations/Opinions of a Johnny Come Lately

Why not smoothen/regularize/parametrize?

• If scale too small it makes no algorithmic difference.
• If scale too large the problem is significantly changed.
• If scale variable we have one more loose parameter.

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 3 / 43

Observations/Opinions of a Johnny Come Lately

Why not smoothen/regularize/parametrize?

• If scale too small it makes no algorithmic difference.
• If scale too large the problem is significantly changed.
• If scale variable we have one more loose parameter.
• Resulting solution paths may have sharp turns (LOP).

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 3 / 43

Observations/Opinions of a Johnny Come Lately

Why not smoothen/regularize/parametrize?

• If scale too small it makes no algorithmic difference.
• If scale too large the problem is significantly changed.
• If scale variable we have one more loose parameter.
• Resulting solution paths may have sharp turns (LOP).
• There may be spurious solutions (Bernardo et al).

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 3 / 43

Observations/Opinions of a Johnny Come Lately

Why not smoothen/regularize/parametrize?

• If scale too small it makes no algorithmic difference.
• If scale too large the problem is significantly changed.
• If scale variable we have one more loose parameter.
• Resulting solution paths may have sharp turns (LOP).
• There may be spurious solutions (Bernardo et al).

Lesson/Moral: Let’s face the combinatorial music! (Reflected in the piecewise linearization)

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 3 / 43

Observations/Opinions of a Johnny Come Lately

My issues with generalized differentiation

• For composite nonsmoothness differentiation rules are only inclusions

so analysis cannot be turned into algebra in an AD fashion.

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 4 / 43

Observations/Opinions of a Johnny Come Lately

My issues with generalized differentiation

• For composite nonsmoothness differentiation rules are only inclusions

so analysis cannot be turned into algebra in an AD fashion.

• The inclusions point the right way for propagating semi-smoothness

but the wrong way for calculating generalized derivative elements.

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 4 / 43

Observations/Opinions of a Johnny Come Lately

My issues with generalized differentiation

• For composite nonsmoothness differentiation rules are only inclusions

so analysis cannot be turned into algebra in an AD fashion.

• The inclusions point the right way for propagating semi-smoothness

but the wrong way for calculating generalized derivative elements.

• Semi-smooth Newton extremely local, bundle methods too heuristic.

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 4 / 43

Observations/Opinions of a Johnny Come Lately

My issues with generalized differentiation

• For composite nonsmoothness differentiation rules are only inclusions

so analysis cannot be turned into algebra in an AD fashion.

• The inclusions point the right way for propagating semi-smoothness

but the wrong way for calculating generalized derivative elements.

• Semi-smooth Newton extremely local, bundle methods too heuristic.
• Generalized derivatives are fickle since outer semi-continuity of

multi-functions does not mean stability in any numerical sense.

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 4 / 43

Observations/Opinions of a Johnny Come Lately

My issues with generalized differentiation

• For composite nonsmoothness differentiation rules are only inclusions

so analysis cannot be turned into algebra in an AD fashion.

• The inclusions point the right way for propagating semi-smoothness

but the wrong way for calculating generalized derivative elements.

• Semi-smooth Newton extremely local, bundle methods too heuristic.
• Generalized derivatives are fickle since outer semi-continuity of

multi-functions does not mean stability in any numerical sense.

• Rademacher says F ∈ C 0,1(Rn) ≡ W 1,∞(Rn), hence

generalized derivatives are almost everywhere normal derivatives.

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 4 / 43

Observations/Opinions of a Johnny Come Lately

Lurking in the background: Prof. Moriarty

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 5 / 43

Piecewise Linearization/Differentiation

Basic idea of tangent linearization:

˚ x F1 F2 F = max(F1, F2) ˚ F1 + F ′

1(˚

x ) ∆ x ˚ F

2

+ F

′ 2

( ˚ x ) ∆ x ˚ F + ∆ F ( ˚ x ; ∆ x )

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 6 / 43

Piecewise Linearization/Differentiation

abs covers min, max and table look-ups

Provided u and w are both finite one has max(u, w) =

1 2 [u + w + abs(u − w)]

min(u, w) =

1 2 [u + w − abs(u − w)]

data (xi, yi) for 0 ≤ i ≤ n are piecewise linearly interpolated by the formula y =

1 2 [y0 + s1 abs(x − x0) + yn + sn abs(x − xn)

+

n−1

• i=1

(si+1 − si) abs(x − xi)] whose ??? where si = (yi+1 − yi)/(xi+1 − xi) represent the slopes.

• Every continuous PL function can be expressed as composition of

affine functions and several abs(). That representation is not unique.

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 7 / 43

Piecewise Linearization/Differentiation

Piecewise Linearization

We wish to determine for base point x and increment ∆x ∆y ≡ ∆F(x; ∆x) = F(x + ∆x) − F(x) + O(∆x2) This can be done by propagating increments according to Smooth elementals ∆vi = ∆vj ± ∆vk for vi = vj ± vk ∆vi = vj ∗ ∆vk + ∆vj ∗ vk for vi = vj ∗ vk ∆vi = cij ∆vj with cij ≡ ϕ′

i(vj)

for vi = ϕi(vj) ≡ abs() Lipschitz Elementals ∆vi = abs(vj + ∆vj) − abs(vj) when vi = abs(vj) . and correspondingly for max() und min().

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 8 / 43

Piecewise Linearization/Differentiation

Continuous Piecewise Differentiation Rules

Linearity and Product Rule F, G : D ⊂ Rn → Rm, α, β ∈ R = ⇒ ∆[αF + βG](x; ∆x) = α∆F(x, ∆x) + β∆G(x, ∆x) ∆[F ⊤G](x; ∆x) = G(x)⊤∆F(x, ∆x) + F(x)⊤∆G(x, ∆x)

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 9 / 43

Piecewise Linearization/Differentiation

Continuous Piecewise Differentiation Rules

Linearity and Product Rule F, G : D ⊂ Rn → Rm, α, β ∈ R = ⇒ ∆[αF + βG](x; ∆x) = α∆F(x, ∆x) + β∆G(x, ∆x) ∆[F ⊤G](x; ∆x) = G(x)⊤∆F(x, ∆x) + F(x)⊤∆G(x, ∆x) Chain Rule F : D ⊂ Rn → Rm and G : E ⊂ Rm → Rp with F(D) ⊂ E = ⇒ ∆[G ◦ F](x; ∆x) = ∆G(F(x); ∆F(x, ∆x))

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 9 / 43

Piecewise Linearization/Differentiation

Second order error and Lipschitz continuity

Proposition Suppose F is composite Lipschitz on some open neighborhood D of a closed convex domain K ⊂ Rn. Then there exists a constant γ such that for all pairs x, x + ∆x ∈ K F(x + ∆x) − F(x) − ∆F(x; ∆x) ≤ γ∆x2

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 10 / 43

Piecewise Linearization/Differentiation

Second order error and Lipschitz continuity

Proposition Suppose F is composite Lipschitz on some open neighborhood D of a closed convex domain K ⊂ Rn. Then there exists a constant γ such that for all pairs x, x + ∆x ∈ K F(x + ∆x) − F(x) − ∆F(x; ∆x) ≤ γ∆x2 Moreover, for any pair ˜ x, x ∈ K, ∆x ∈ Rn, and a constant ˜ γ ∆F(˜ x; ∆x) − ∆F(x; ∆x) /(1 + ∆x) ≤ ˜ γ˜ x − x

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 10 / 43

Piecewise Linearization/Differentiation

Second order error and Lipschitz continuity

Proposition Suppose F is composite Lipschitz on some open neighborhood D of a closed convex domain K ⊂ Rn. Then there exists a constant γ such that for all pairs x, x + ∆x ∈ K F(x + ∆x) − F(x) − ∆F(x; ∆x) ≤ γ∆x2 Moreover, for any pair ˜ x, x ∈ K, ∆x ∈ Rn, and a constant ˜ γ ∆F(˜ x; ∆x) − ∆F(x; ∆x) /(1 + ∆x) ≤ ˜ γ˜ x − x Finally there is a continuous radius ρ(x) such that ∆F(x; ∆x) = F ′(x; ∆x) if ∆x < ρ(x) Locally we reduce to the homogeneous piecewise linear F ′(x; ∆x).

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 10 / 43

Piecewise Linearization/Differentiation

Differentiation Concepts on Euclidean Spaces

Function Space: Diff.Op.: Model Space: Discrepancy: Smooth = S ∂|˚

x

− → Lip L = linear uniform ⊂ ⊂ CompPS = CPS ∆|˚

x

− → Lip PL = Piecewise L uniform ⊂ − → ∂B

• ˚

x

LipschitzPS = LPS ∂B

• ˚

x

− → ??? PLh = homog. PL nonuniform ⊂ PiecewiseS = DCPS ∆|˚

x

− → ??? DPL = discont. PL nonuniform

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 11 / 43

Piecewise Linearization/Differentiation

Piecewise Linearization of Discontinuous f

where (at the origin)

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 12 / 43

Representation of PL functions in abs-normal form

Polyhedral Decomposition

x 0.0 0.5 1.0 1.5 2.0 y 0.0 0.5 1.0 1.5 2.0 2.5 3.0 F ( x , y ) 2 4 6 8 10 12 14 16 18 Polyhedron

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 13 / 43

Representation of PL functions in abs-normal form

A simple R2 → R2 example:

f1 = x1 + |x1 − x2| + |x1 − |x2||, f2 = x2 The switching variables zi are the arguments of the abs-functions: z1 = x1 − x2, z2 = x2, z3 = x1 − |z2| ⇒ f1 = x1 + |z1| + |z3| In Abs-normal form:       z1 z2 z3 y1 y2       =      

Z

1 −1

L

• 1

1 −1 1 1 1

J

1

• Y

            x1 x2 |z1| |z2| |z3|      

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 14 / 43

Representation of PL functions in abs-normal form

Graph reprsentation and switching depth

Computational graph

x1 z1 z3 z2 y1 x2 y2

Proposition (Griewank et al 2014) Switching depth ν ≡ length of directed path is bounded by ¯ ν(n) = 2 n − 1 Example: ordervaluek(z1, . . . , zm) ≡ k−th largest zi (N. Kreijic)

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 15 / 43

Representation of PL functions in abs-normal form

Abs-normal form for y = F(x) : Rn → Rm

z y

• =

c b

• +

Z L J Y x |z|

• J ∈ Rm×n represents the smooth part of the function
• L ∈ Rs×s is strictly lower triangular to yields z = z(x) uniquely.
• Crossterms Z ∈ Rs×n and Y ∈ Rm×s link x to z and z to y.

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 16 / 43

Representation of PL functions in abs-normal form

Abs-normal form for y = F(x) : Rn → Rm

z y

• =

c b

• +

Z L J Y x |z|

• J ∈ Rm×n represents the smooth part of the function
• L ∈ Rs×s is strictly lower triangular to yields z = z(x) uniquely.
• Crossterms Z ∈ Rs×n and Y ∈ Rm×s link x to z and z to y.
• The switching depth ν is the structural nilpotency degree of L.
• Abs-normal form is nonredundant and stable w.r.t. perturbations.
• ADOL-C can calculate [b, c, Z, L, J, Y ] after slight modification.

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 16 / 43

Representation of PL functions in abs-normal form

Abs-normal form for y = F(x) : Rn → Rm

z y

• =

c b

• +

Z L J Y x |z|

• J ∈ Rm×n represents the smooth part of the function
• L ∈ Rs×s is strictly lower triangular to yields z = z(x) uniquely.
• Crossterms Z ∈ Rs×n and Y ∈ Rm×s link x to z and z to y.
• The switching depth ν is the structural nilpotency degree of L.
• Abs-normal form is nonredundant and stable w.r.t. perturbations.
• ADOL-C can calculate [b, c, Z, L, J, Y ] after slight modification.
• The sign vector σ ≡ sign(z(x)) determines the control flow.
• Σ ≡ diag(σ) =

⇒ |z| = Σz is componentwise modulus.

• Relatively open Pσ = {x : σ(x) = σ} form polyhedral skeleton.

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 16 / 43

Representation of PL functions in abs-normal form

Some other References and Traditions

Stern, Thomas Edwin (1956) Piecewise-linear network theory

• W. M. G. van Bokhoven (1981)

Piecewise-linear modelling and analysis

• J. Rohn (2004)

A Theorem of the Alternatives for the equation Ax + B|x| = b

• O. Mangasarian, and R. Meyer( 2006)

Absolute value equations

• L. Brugnano, and V. Casulli (2008)

Iterative solution of piecewise linear systems

• A. Martin (2009)

Mixed Integer Nonlinear Programming

• A. Hadjidimos, M. Lapidakis, M. and M. Tzoumas (2012)

On Iterative Solution for Linear Complementarity Problem with an H+-Matrix

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 17 / 43

Computation of conical Jacobians and gradients

Limiting Jacobians/Gradients of PL function

Limiting (=Bouligand) Jacobian: ∂LF(x) ≡ {Jσ|x ∈ ¯ Pσ, with Pσ open} where Jσ = J + Y Σ(I − LΣ)−1Z with (I − LΣ)−1 = I + LΣ + (LΣ)2 + · · · + (LΣ)(ν−1)

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 18 / 43

Computation of conical Jacobians and gradients

Limiting Jacobians/Gradients of PL function

Limiting (=Bouligand) Jacobian: ∂LF(x) ≡ {Jσ|x ∈ ¯ Pσ, with Pσ open} where Jσ = J + Y Σ(I − LΣ)−1Z with (I − LΣ)−1 = I + LΣ + (LΣ)2 + · · · + (LΣ)(ν−1) Polynomial escape in direction d1 x(t) = x +

n

• i=1

tidi ∈ Pσ open when det(d1 . . . dn) = 0 and 0 < t ≈ 0

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 18 / 43

Computation of conical Jacobians and gradients

Limiting Jacobians/Gradients of PL function

Limiting (=Bouligand) Jacobian: ∂LF(x) ≡ {Jσ|x ∈ ¯ Pσ, with Pσ open} where Jσ = J + Y Σ(I − LΣ)−1Z with (I − LΣ)−1 = I + LΣ + (LΣ)2 + · · · + (LΣ)(ν−1) Polynomial escape in direction d1 x(t) = x +

n

• i=1

tidi ∈ Pσ open when det(d1 . . . dn) = 0 and 0 < t ≈ 0 Complexity range (also) utilizing reverse mode 3 min(m, n) ≤ OPS(Jσ)/OPS(F) ≤ 3 n

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 18 / 43

Computation of conical Jacobians and gradients

Conical Jacobians of PS function

Proposition: Khan & Barton and A. G. 2013 ∂KF(x) ≡ ∂L

∆x∆F(x; ∆x)

• ∆x=0 ⊂ ∂LF(x)

contains exactly those Jacobians ∂Fσ(x) for which the tangent cone Tσ ≡ Tx{z ∈ D : Fσ(z) = F(z)} has a nonempty interior. (i.e. Fσ and ∂Fσ are conically active)

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 19 / 43

Computation of conical Jacobians and gradients

Conical Jacobians of PS function

Proposition: Khan & Barton and A. G. 2013 ∂KF(x) ≡ ∂L

∆x∆F(x; ∆x)

• ∆x=0 ⊂ ∂LF(x)

contains exactly those Jacobians ∂Fσ(x) for which the tangent cone Tσ ≡ Tx{z ∈ D : Fσ(z) = F(z)} has a nonempty interior. (i.e. Fσ and ∂Fσ are conically active) Kummer example for ∂Lf = ∂Kf ≡ ∂L∆f (x; ·)

−1 −0.5 0.5 1 −1 −0.5 0.5 1 0.2 0.4 0.6 0.8 1

f (x, y) = (y2 − x+)+ with z+ ≡ max(0, z) = ⇒ {0} = ∂Kf (0) ∋ (−1, 0) ∈ ∂Lf (0) ⊂ ∂Cf (0)

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 19 / 43

Solving PL systems of Equation and LCPs

Related perturbed semismooth system

F(x, y) = x/2 + (y 2 − x+)+ y

• =

δ ε

• has for x ≤ 0, 0 < x < y 2, and y 2 < x, respectively

F ′(x, y) = 1/2 2y 1

• ,

−1/2 2y 1

• and

1/2 1

• Figure: Cyclic behavior of Newton = Sucessive Linearization

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 20 / 43

Solving PL systems of Equation and LCPs

Reformulations in the equation case m = n

Original Equation 0 = F(x) = b + J x + Y z(x) with z(x) = (I − LΣ)−1(c + Z x)

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 21 / 43

Solving PL systems of Equation and LCPs

Reformulations in the equation case m = n

Original Equation 0 = F(x) = b + J x + Y z(x) with z(x) = (I − LΣ)−1(c + Z x) Schur complement and complementary system det(J) = 0 ( achievable using v ≡ |v + |v|| − |v|) ensures existence S ≡ L − ZJ−1Y ∈ Rs×s = ⇒ z = S|z| + ˆ c with ˆ c = c − Z J−1b

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 21 / 43

Solving PL systems of Equation and LCPs

Reformulations in the equation case m = n

Original Equation 0 = F(x) = b + J x + Y z(x) with z(x) = (I − LΣ)−1(c + Z x) Schur complement and complementary system det(J) = 0 ( achievable using v ≡ |v + |v|| − |v|) ensures existence S ≡ L − ZJ−1Y ∈ Rs×s = ⇒ z = S|z| + ˆ c with ˆ c = c − Z J−1b Corresponding Linear Complementarity Problem z = u − w, 0 ≤ u ⊥ w ≥ 0 = ⇒ u ⊥ Mu + q ≥ 0 with M = (I − S)−1(I + S)

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 21 / 43

Solving PL systems of Equation and LCPs

Conditions for solvability and convergence

General Implication chain for PL functions

Bijectivity ⇔ Injectivity ⇒ Openness ⇔ Coherent Orientation ⇒ Surjectivity

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 22 / 43

Solving PL systems of Equation and LCPs

Conditions for solvability and convergence

General Implication chain for PL functions

Bijectivity ⇔ Injectivity ⇒ Openness ⇔ Coherent Orientation ⇒ Surjectivity

Griewank et al 2013: For equation 0 = F(x) = min(x, Mx + q) and other simply switched system coherent orientation (i.e. M is P-Matrix) implies already bijectivity.

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 22 / 43

Solving PL systems of Equation and LCPs

Conditions for solvability and convergence

General Implication chain for PL functions

Bijectivity ⇔ Injectivity ⇒ Openness ⇔ Coherent Orientation ⇒ Surjectivity

Griewank et al 2013: For equation 0 = F(x) = min(x, Mx + q) and other simply switched system coherent orientation (i.e. M is P-Matrix) implies already bijectivity. Sherman-Morrison-Woodbury yields det(Jσ) = det(J) det(I − S Σ) Hence det(I − SΣ) > 0 for all Σ ⇒ Coherent Orientation

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 22 / 43

Solving PL systems of Equation and LCPs

Equivalent and related conditions

Rump showed With sign real spectral radius ρs

0(S) = max(|λ|) over R ∋ λ ∈ spect(S)

det(I − S Σ) > 0 ⇔ ρs

0(S) < 1

⇔ (S − I)−1(S + I) is P also equivalent to non-expansiveness (Rohn). x ≥ 0 = ⇒ |S x| ≥ |x| componentwise Implication chain Difficulty: Test for above property is NP-Hard. ρ(|S|) < 1 ⇒ D−1SDp < 1 ⇒ ρs

0(S) < 1

(Absolute contractivity) (Smooth dominance) (Coherent orientation.)

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 23 / 43

Solving PL systems of Equation and LCPs

Solvers for PL systems of equations in

• riginal abs-normal or complementary form.

Method Convergence condition Rate Effort Generalized Newton on OPL 2 ˆ ρ < (1−Lp−ˆ ρ/2)2 finite I −SΣ, J Generalized Newton on CPL Sp < 1/3 finite I −SΣ Signed Gauss on CPL ρ(|S|) ≤ 1/2 finite I −SΣ once Block Seidel on CPL S − Lp + Lp < 1 linear I −LΣ, J Modulus Iteration on CPL Sp < 1 linear J Piecewise Newton on OPL coherent orient. of F finite I −SΣ, J Piecewise Newton on CPL ρs

0(S) < 1

finite I −SΣ

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 24 / 43

Solving PL systems of Equation and LCPs

Numerical results on Murty’s LCP example

n Anzahl der Iterationen PLN (1. Startvektor) PLN (2. Startvektor) Lemke 2 2 2 4 = 22 4 4 6 16 = 24 6 10 14 64 = 26 8 24 32 256 = 28 10 54 66 1024 = 210 12 116 136 4096 = 212 14 246 276 16384 = 214 16 512 558 216 18 1048 1110 218 20 2126 2220 220

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 25 / 43

Solving PL systems of Equation and LCPs

Piecewise linear Newton on Rosette Example

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 26 / 43

(Un)constrained optimization by successive PL

Optimization with quadratic overestimation

Under our assumptions on compact domains ˆ q(x, ∆x) ≡ |f (x + ∆x) − f (x) − ∆f (x; ∆x)| ∆x2 ≤ ¯ q(∆x) Consequence: Bundle type iteration ∆x ≡ argmin

˜ s

(∆f (x;˜ s) + q˜ s2) x += ∆x if f (x + ∆x) < f (x) q+ = max(q, ˆ q(x, ∆x)) is guaranteed to converge from within bounded level set.

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 27 / 43

(Un)constrained optimization by successive PL

Application to Rosenbrock ´ a la Nesterov

f (x1, x2) = 1

4(x1 − 1)2 +

• x2 − 2 x2

1 + 1

• .

yields piecewise linearization f (x1, x2)+∆f (x1, x2; ∆x1, ∆x2) =

1 4(x1 − 1)2 + 1 2(x1 − 1)∆x1 +

• x2 + ∆x2 − 2 x2

1 − 4x1∆x1 + 1

• .

x2 x1 −2 −1 1 2 3 −2 −1 1 2 3 contours starting point 1 starting point 2 starting point 3

100 200 300 400 500 10

−6

10

−4

10

−2

10 10

2

iterations function value starting point 1 starting point 2 starting point 3

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 28 / 43

(Un)constrained optimization by successive PL

Local = Inner Problem

min

s∈Rn ∆f (x; s) + q 2s2

• At least, global minimization is NP-hard ( ← SAT3)
• Bad News going back to Hirriart-Urruty & Lemarechal:

Steepest descent with exact line search may fail on convex PL f .

• Challenge is to avoid Zenon effect = Zigzagging

−100 −50 50 −20 −15 −10 −5 5 10 15 20

x1 x2 Nondifferentiable points of f f0(x) f2(x) f−2(x) f1(x) f−1(x) x0=(9,−3) Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 29 / 43

(Un)constrained optimization by successive PL

Good News by H. U. & L and Griewank et al:

True steepest descent trajectory x(t) defined by: −d x(t) d t+ = −d(x) ≡ short(∂f (x)) ≡ argmin{g : g ∈ ∂f (x)} is in convex case unique solution of differential inclusion ˙ x ∈ −∂f (x), which has stationary cluster points or limit x∗ in initial level set. Can be realized using abs-normal form an and Zenon effect excluded.

−100 −50 50 −20 −15 −10 −5 5 10 15 20

x1 x2 Nondifferentiable points of f f0(x) f2(x) f−2(x) f1(x) f−1(x) x0=(9,−3) Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 30 / 43

(Un)constrained optimization by successive PL

Good News by H. U. & L and Griewank et al:

True steepest descent trajectory x(t) defined by: −d x(t) d t+ = −d(x) ≡ short(∂f (x)) ≡ argmin{g : g ∈ ∂f (x)} is in convex case unique solution of differential inclusion ˙ x ∈ −∂f (x), which has stationary cluster points or limit x∗ in initial level set. Can be realized using abs-normal form an and Zenon effect excluded.

−100 −50 50 −20 −15 −10 −5 5 10 15 20

x1 x2 Nondifferentiable points of f f0(x) f2(x) f−2(x) f1(x) f−1(x) x0=(9,−3) Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 30 / 43

(Un)constrained optimization by successive PL

Good News by H. U. & L and Griewank et al:

True steepest descent trajectory x(t) defined by: −d x(t) d t+ = −d(x) ≡ short(∂f (x)) ≡ argmin{g : g ∈ ∂f (x)} is in convex case unique solution of differential inclusion ˙ x ∈ −∂f (x), which has stationary cluster points or limit x∗ in initial level set. Can be realized using abs-normal form an and Zenon effect excluded.

−100 −50 50 −20 −15 −10 −5 5 10 15 20

x1 x2 Nondifferentiable points of f f0(x) f2(x) f−2(x) f1(x) f−1(x) x0=(9,−3) Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 30 / 43

(Un)constrained optimization by successive PL

Good News by H. U. & L and Griewank et al:

True steepest descent trajectory x(t) defined by: −d x(t) d t+ = −d(x) ≡ short(∂f (x)) ≡ argmin{g : g ∈ ∂f (x)} is in convex case unique solution of differential inclusion ˙ x ∈ −∂f (x), which has stationary cluster points or limit x∗ in initial level set. Can be realized using abs-normal form an and Zenon effect excluded.

−100 −50 50 −20 −15 −10 −5 5 10 15 20

x1 x2 Nondifferentiable points of f f0(x) f2(x) f−2(x) f1(x) f−1(x) x0=(9,−3) x*=(−50,0) Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 30 / 43

(Un)constrained optimization by successive PL −200 −150 −100 −50 50 100 −50 −40 −30 −20 −10 10 20 30 40 50

x1 x2 Proximal Bundle Method f0(x) f2(x) f−2(x) f1(x) f−1(x) x0=(9,−3) x*=(−157.24,−19.62)

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 31 / 43

Integration of Lipschitzian dynamics

ODE integration with Lipschitzian RHS

Possibly after space discretization of PDE: ˙ x ≡ d dt x(t) = F(x(t)) with F ∈ C0,1 = W 1,∞

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 32 / 43

Integration of Lipschitzian dynamics

ODE integration with Lipschitzian RHS

Possibly after space discretization of PDE: ˙ x ≡ d dt x(t) = F(x(t)) with F ∈ C0,1 = W 1,∞ Generalized midpoint rule With ˇ x current point, ˆ x next point, ˚ x = (ˇ x + ˆ x)/2 and time step h ˆ x − ˇ x = h 1/2

−1/2

[F(˚ x) + ∆F (˚ x; (ˆ x − ˇ x) t)] dt yields local third order truncation and globally second order convergence.

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 32 / 43

Integration of Lipschitzian dynamics

ODE integration with Lipschitzian RHS

Possibly after space discretization of PDE: ˙ x ≡ d dt x(t) = F(x(t)) with F ∈ C0,1 = W 1,∞ Generalized midpoint rule With ˇ x current point, ˆ x next point, ˚ x = (ˇ x + ˆ x)/2 and time step h ˆ x − ˇ x = h 1/2

−1/2

[F(˚ x) + ∆F (˚ x; (ˆ x − ˇ x) t)] dt yields local third order truncation and globally second order convergence. Properties in special cases When F is PL GMP coincides with Average Vector Field Method (Quispel). Thus exact energy preservation if F = J ∇f Hamiltonian. Generally GMP is in contrast to IMP not (nonsmooth) symplectic.

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 32 / 43

Integration of Lipschitzian dynamics

Experiments Chua circuit

Problem Definition

F(x) =

  ˙

x

˙

y

˙

z

  =   α(y − x − f(x))

x − y + z

−βy  

f(x) = m1x + 1 2(m0 − m1)(|x + 1|−|x − 1|) x,y are the voltages across C1 and C2 z is the intensity of the electrical current at I f(x) is the electrical response of the resistor constants are α = 15.6,β = 28,m0 = −1.143,m1 = −0.714 .

Figure: Chua circuit

taken from

http://www.chuacircuits.com/

Paul Boeck Lipschitzean ODEs March 20, 2013 17 / 22

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 33 / 43

Integration of Lipschitzian dynamics

Experiments Chua circuit

Chua circuit

−4 −2 2 4 −0.4 −0.2 0.2 0.4 −4 −2 2 4 x y z

Figure: Chua Circuit - The Double Scroll

Paul Boeck Lipschitzean ODEs March 20, 2013 16 / 22

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 34 / 43

Integration of Lipschitzian dynamics

Experiments Chua circuit

Convergence

10

2

10

3

10

4

10

5

10

−7

10

−6

10

−5

10

−4

10

−3

10

−2

10

−1

10 10

1

Number of steps Error in x New Midpoint Old Midpoint

Figure: Error compared to fine grid solution

Paul Boeck Lipschitzean ODEs March 20, 2013 20 / 22

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 35 / 43

Integration of Lipschitzian dynamics

First level Richardson Extrapolation yields

103 104 10−10 10−9 10−8 10−7 10−6 10−5 10−4 10−3 10−2 4 2 degrees of freedom N error over all components Generalized Midpoint Implicit Midpoint Explicit Midpoint

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 36 / 43

Integration of Lipschitzian dynamics

Electrical circuit with diode

V (t) g(x) C L I(t) = z(t)

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 37 / 43

Integration of Lipschitzian dynamics

0.5 1 1.5 ·10−8 1 2 3 ·10−4 I(t) 0.5 1 1.5 ·10−8 0.5 1 1.5 ·10−13 Q(t)

Figure: Solution of the ODE System

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 38 / 43

Integration of Lipschitzian dynamics

Log-log plot of convergence

102 103 104 105 10−10 10−9 10−8 10−7 10−6 10−5 10−4 10−3 Number of timesteps n Error in x

• Impl. Midpoint

Trapezoidal

• Gen. Midpoint
• Gen. Trapezoidal

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 39 / 43

Integration of Lipschitzian dynamics

Richardson/Romberg Extrapolation

104 105 106 107 Number of steps 10-13 10-12 10-11 10-10 10-9 Error

Convergence

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 40 / 43

Integration of Lipschitzian dynamics

Convergence Orders

Midpoint/Trapezoidal Classical Generalized Extrapolated General Position O(h) O(h2) O(h2) Transversal Position O(h2) ch2 + O(h3) O(h3) PL+smooth forcing O(h2) ch2 + O(h4) O(h4) Conjecture: On discontinuous right hand sides all orders reduced by 1. Challenge: Algebraic inclusion solving In discontinuous case PL based discretization requires solution of F(x) ∋ 0. Special case: piecewise constant F = ∇f for continuous PL objective f .

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 41 / 43

Integration of Lipschitzian dynamics

Summary and Conclusion

• Differentiation of composite piecewise smooth functions yields PL

function in abs-normal form.

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 42 / 43

Integration of Lipschitzian dynamics

Summary and Conclusion

• Differentiation of composite piecewise smooth functions yields PL

function in abs-normal form.

• Local PL model has uniform second order fit and varies Lipschitz

continuously with reference point.

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 42 / 43

Integration of Lipschitzian dynamics

Summary and Conclusion

• Differentiation of composite piecewise smooth functions yields PL

function in abs-normal form.

• Local PL model has uniform second order fit and varies Lipschitz

continuously with reference point.

• Limiting Jacobians/gradients of PL model can be computed and are

conical derivatives of CPS function.

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 42 / 43

Integration of Lipschitzian dynamics

Summary and Conclusion

• Differentiation of composite piecewise smooth functions yields PL

function in abs-normal form.

• Local PL model has uniform second order fit and varies Lipschitz

continuously with reference point.

• Limiting Jacobians/gradients of PL model can be computed and are

conical derivatives of CPS function.

• When m = n the PL model system can be solved directly or after

rewriting as complementary problem. In the transformation openness = coherent orientation = path continuation may be lost.

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 42 / 43

Integration of Lipschitzian dynamics

Summary and Conclusion

• Differentiation of composite piecewise smooth functions yields PL

function in abs-normal form.

• Local PL model has uniform second order fit and varies Lipschitz

continuously with reference point.

• Limiting Jacobians/gradients of PL model can be computed and are

conical derivatives of CPS function.

• When m = n the PL model system can be solved directly or after

rewriting as complementary problem. In the transformation openness = coherent orientation = path continuation may be lost.

• Abs-normal form provides grey box information for systematic
• ptimization of PL function by bundle type method.

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 42 / 43

Integration of Lipschitzian dynamics

Summary and Conclusion

• Differentiation of composite piecewise smooth functions yields PL

function in abs-normal form.

• Local PL model has uniform second order fit and varies Lipschitz

continuously with reference point.

• Limiting Jacobians/gradients of PL model can be computed and are

conical derivatives of CPS function.

• When m = n the PL model system can be solved directly or after

rewriting as complementary problem. In the transformation openness = coherent orientation = path continuation may be lost.

• Abs-normal form provides grey box information for systematic
• ptimization of PL function by bundle type method.
• Lipschitzian dynamical systems can be integrated with generalized

and extrapolated midpoint/trapezoidal rule of global order 2, 3 or 4.

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 42 / 43

Integration of Lipschitzian dynamics

Summary and Conclusion

• Differentiation of composite piecewise smooth functions yields PL

function in abs-normal form.

• Local PL model has uniform second order fit and varies Lipschitz

continuously with reference point.

• Limiting Jacobians/gradients of PL model can be computed and are

conical derivatives of CPS function.

• When m = n the PL model system can be solved directly or after

rewriting as complementary problem. In the transformation openness = coherent orientation = path continuation may be lost.

• Abs-normal form provides grey box information for systematic
• ptimization of PL function by bundle type method.
• Lipschitzian dynamical systems can be integrated with generalized

and extrapolated midpoint/trapezoidal rule of global order 2, 3 or 4.

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 42 / 43

Integration of Lipschitzian dynamics

References

• A.G.: On stable piecewise linearization and generalized algorithmic
• differentiation. Optimization Methods and Software, (2013).
• P. Boeck, B. Gompil, A. Griewank, R. Hasenfelder, N. Strogies:

Experiments with Generalized Midpoint and Trapezoidal Rules on two Nonsmooth ODEs, Mongolian Journal of Mathematics (2014).

• A.G., M. Radons, T. Streubel, J.U. Bernd:

On piecewise linear equations in abs-normal form. Revised submission to Linear Algebra and Applications (2014).

• A.G., A. Walther, S. Fiege, T. Bosse: On Lipschitz optimization

based on gray-box piecewise linearization submitted to Mathematical Programming, (2014)

• A.G. et al: Plan-C, Piecewise Linear functions in Abs Normal form,

software under development.

Andreas Griewank et al (HUB) Nonsmooth numerics via PL August 8, 2014 43 / 43