Generalized Derivatives Automatic Evaluation & Implications for - - PowerPoint PPT Presentation

Generalized Derivatives

Automatic Evaluation & Implications for Algorithms Paul I. Barton, Kamil A. Khan & Harry A. J. Watson

Process Systems Engineering Laboratory Massachusetts Institute of Technology

2

Nonsmooth Equation Solving

◆ Semismooth Newton method: ◆ Linear programming (LP) Newton method: ◆ some element of a generalized derivative

min

γ ,x γ

s.t. f(xk ) + G(xk )(x − xk )

∞ ≤ γ f(xk ) ∞ 2

(x − xk )

∞ ≤ γ f(xk ) ∞

x ∈X

G(xk )(x − xk ) = −f(xk )

Kojima & Shindo (1986), Qi & Sun (1993), Facchinei, Fischer & Herrich (2014).

G(xk )

Polyhedral set

3

◆ Suppose locally Lipschitz => differentiable on a set S ◆ B-subdifferential: ◆ Clarke Jacobian: ◆ Useful properties of :

Ø Nonempty, convex, and compact Ø Satisfies mean-value theorem, implicit/inverse function theorems Ø Reduces to subdifferential/derivative when is convex/strictly differentiable

Generalized Derivatives

∂Bf(x):={H :H = lim

i→∞Jf(x(i)), x = lim i→∞x(i), x(i) ∈S}

∂f(x):= conv∂Bf(x)

f (x) = x ∂f (x) = {1} ∂f (x) = {−1} x ∂B f (x) ={−1,1}, ∂ f (x) = [−1,1]

∂f(x)

f

f

Clarke (1973).

4

Convergence Properties

◆ Suppose generalized derivative contains no

singular matrices at the solution

◆ Semismooth Newton method: Ø local Q-superlinear convergence if Ø local Q-quadratic convergence if strongly semismooth ◆ Semismooth Newton & LP-Newton methods for

PC1 or strongly semismooth functions:

Ø local Q-quadratic convergence if ◆ Automatic/Algorithmic Differentiation (AD) Ø Automatic methods for computing derivatives in complex settings Ø Automatic method for computing elements of generalized derivatives? Ø Computationally relevant generalized derivatives

G(xk ) ∈∂Bf(xk ) G(xk ) ∈∂f(xk )

5

All generalized derivatives are equal… But, some are more equal than others.

6

Obstacles to Automatic Gen. Derivative Evaluation 1

◆ Automatically evaluating Clarke Jacobian

elements is difficult

◆ Lack of sharp calculus rules: g(x) = max{0, x} 0 ∈ ∂g(0) =[0,1]

x x x

0 ∈ ∂h(0) =[0,1] (0 + 0) ∉∂ f (0) ={1} h(x) = min{0, x} f (x) = g(x)+ h(x) ∂ f (0) ⊂ ∂g(0) + ∂h(0)

7

◆ Directional derivative: ◆ Sharp chain rule for locally Lipschitz functions: ◆ AD gives the directional derivative ◆ PC1 functions: finite collection of C1 functions for which ◆ 2-norm not PC1

Directional Derivatives & PC1 Functions

f '(x;d) = lim

t→0+

f(x + td) − f(x) t

[f !g]'(x;d) = f '(g(x);g'(x;d))

Griewank (1994), Scholtes (2012).

f(y) ∈ φ(y):φ ∈Ff (x)

{ }, ∀y ∈N(x)

8

Obstacles 2

◆ PC1 functions have piecewise linear directional

derivative

d1 d2

′ f (x;d) = B(1) d ′ f (x;d) = B(2) d ′ f (x;d) = B(3) d

9

Obstacles 2

◆ PC1 functions have piecewise linear directional

derivative

◆ Directional derivatives in the coordinate

directions do not necessarily give B- subdifferential elements

◆ Also defeats finite differences

d1 d2

′ f (x;d) = B(1) d ′ f (x;d) = B(2) d ′ f (x;d) = B(3) d

10

Obstacles 3

◆ may be a strict subset of

∂f(x)

∂ fi(x)

i=1 m

∏

f :(x1, x2) ! x1+ | x2 | x1− | x2 | ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ ∂f(0) = 1 1 ⎡ ⎣ ⎢ 2s −1 1− 2s ⎤ ⎦ ⎥ : s ∈ 0,1 ⎡ ⎣ ⎤ ⎦ ⎧ ⎨ ⎩ ⎪ ⎫ ⎬ ⎭ ⎪ ∂ f1(0) × ∂ f2(0) = 1 1 ⎡ ⎣ ⎢ 2s1 −1 2s2 −1 ⎤ ⎦ ⎥ :(s1,s2) ∈ 0,1 ⎡ ⎣ ⎤ ⎦

2

⎧ ⎨ ⎪ ⎩ ⎪ ⎫ ⎬ ⎪ ⎭ ⎪

π 2 ∂f(0) π 2(∂ f1(0) × ∂ f2(0))

11

L-smooth Functions

Nesterov (1987), Khan and Barton (2014), Khan and Barton (2015).

◆ The following functions are L-smooth:

Ø Continuously differentiable functions Ø Convex functions (e.g. abs, 2-norm) Ø PC1 functions Ø Compositions of L-smooth functions: Ø Integrals of L-smooth functions: Ø Solutions of ODEs with L-smooth right-hand sides: where

f : X ∈ Rn → Rm x ! h(g(x)) x ! g(t,x) dt

a b

∫

dx dt (t,c) = g(t,x(t,c)), x(0,c) = c c ! x(b,c),

12 ◆

L-subdifferential:

Ø Contains L-derivatives in directions M:

◆

Useful properties:

Ø L-derivatives classical derivative wherever strictly differentiable Ø L-derivatives elements of Clarke gradient Ø Contains only subgradients when f convex Ø Contained in plenary hull of Clarke Jacobian, and can be used in place of Clarke Jacobian in numerical methods: Ø For PC1 functions, L-derivatives elements of B-subdifferential Ø Satisfies sharp chain rule, expressed naturally using LD-derivatives

Lexicographic Derivatives

{Ad: A ∈∂Lf(x)}⊂{Ad: A ∈∂f(x)} for each d ∈Rn

∂Lf(x) ={JLf(x;M):detM ≠ 0}

JLf(x;M), detM ≠ 0

Nesterov (1987), Khan and Barton (2014), Khan and Barton (2015).

13

◆ Extension of classical directional derivative ◆ LD-derivative: for any ◆ If M is square and nonsingular: ◆ If f is differentiable at x: ◆ Sharp LD-derivative chain rule:

Lexicographic Directional (LD)-Derivatives

M := [m(1) ! m( p)]∈Rn× p,

f '(x;M) = [fx,M

(0) (m(1)) ! fx,M ( p−1)(m( p))]

f '(x;M) = JLf(x;M)M [f !g]'(x;M) = f '(g(x);g'(x;M)) f '(x;M) = Jf(x)M

Khan and Barton (2015).

14

◆ Sharp chain rule immediately implies, given the “seed

directions” M, forward-mode AD can compute:

◆ Need calculus rules for “elementary functions”:

Ø abs, min, max, mid, , etc. Ø algorithm for “elemental PC1 functions” Ø linear programs and lexicographic linear programs parameterized by their RHSs Ø implicit function: is the unique solution N of

Vector Forward AD Mode for LD-derivatives

f '(x;M)

⋅ 2 Khan and Barton (2015), Khan and Barton (2013), Hoeffner et al. (2015).

h(w(z),z) = 0

h' (ˆ y,ˆ z);(N,M)

( ) = 0

w'(ˆ z;M)

15

◆ Inexact Newton method: ◆ Solve iteratively: ◆ But, directional derivative not a linear function of the

directions…

◆ Let , M nonsingular. Then: ◆ But, M not known in advance ◆ Compute columns of one at time

Ø computation of a column affects subsequent columns Ø automatic code can be “locked” to record influence of earlier columns

◆ Local Q-superlinear & Q-quadratic convergence rates

can be achieved

Semismooth Inexact Newton Method

J(x)di, i = 1,2,…

JLf(x;M) Δx = −f(x)

M = d1,d2,… ⎡ ⎣ ⎤ ⎦

f '(x;M) = JLf(x;M)M

f '(x;M)

16

Approximation of LD-derivatives using FDs

◆

LD-derivative:

◆

FD approx. of using p+1 function evaluations: M := [m(1) ! m( p)]∈Rn× p

f '(x;M) = [fx,M

(0) (m(1)) ! fx,M ( p−1)(m( p))]

fx,M

(0) (m(1)) ≈ α −1[f(x +αm(1)) − f(x)] =: Dαm(1)[f](x)

fx,M

(1) (m(2)) ≈ Dαm(2)[fx,M (0) ](m(1)) = Dαm(2)Dαm(1)[f](x)

! fx,M

( p−1)(m( p)) ≈ Dαm( p)[fx,M ( p−2)](m( p−1)) = Dαm( p)"Dαm(2)Dαm(1)[f](x)

x +αm(1) +α 2m(2) x +αm(1) x fx,M

(1) (m(2)) ≈ α −2[f(x +αm(1) +α 2m(2)) − f(x +αm(1))]

f '(x;M)

f '(x;M) = [fx,M

(0) (m(1)) fx,M (1) (m(2))]

fx,M

(0) (m(1)) ≈ α −1[f(x +αm(1)) − f(x)]

17

◆ Cost of AD can be reduced when the Jacobian is sparse

Ø Find structurally orthogonal columns Ø Perform vector forward pass with seed matrix rather than

◆ AD for LD-derivatives à order of the directions matters

Ø Corresponding to M is an uncompressed (permutation) matrix Q:

» M = QD for some matrix D

Ø Procedure: » Identify matrices Q, D, and M » Perform vector forward pass to calculate » Copy entries of into entries of sparse data structure for

◆ Done based on assumption that

» Calculate (i.e. by sparse permutation)

Ø is not true in general a b c d e f g h ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥

Sparse Accumulation for L- derivatives

1 1 1 1 ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ 1 1 1 1 ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ ⎥

M ∈ ϒ

n× p

I ∈ ϒ

n×n

′ f (x;M) ′ f (x;M) ′ f (x;Q) ′ f (x;M) = ′ f (x;Q)D ′ f (x;M) = ′ f (x;Q)D JLf(x;Q) = ′ f (x;Q)Q−1

18

Generalized Derivatives of Algorithms: MHEX model

F

i Ti in − Ti

• ut

( )

i∈H

∑

= f j ti

• ut − ti

in

( )

j∈ C

∑

min

p∈P EBP C p − EBP H p

( ) = 0

UA− ΔQk ΔTLM

k k∈K k≠|K|

∑

= 0 MHEX

F

1,T 1 in

F

1,T 1

• ut

f|C|,t|C|

• ut

f|C|,t|C|

• ut

! F

|H|,T|H|

• ut

F

|H|,T|H|

• ut

! ! ! f1,t1

in

f1,t1

• ut

Watson et al. (2015).

19

◆ ◆ The LD-derivative mapping uniquely

solves the ODE:

◆ Boundary value problem: ◆ Solve with semismooth (inexact) Newton method using

chain rule for LD-derivatives

Ø if in addition f is semismooth, then Q-superlinear convergence rate Ø if it happens to be PC1, then Q-quadratic convergence rate

ODEs and BVPs

dx dt (t,c) = f(t,x(t,c)), x(t0,c) = c

t ! [xt]'(c0;M)

dA dt (t) = [ft]'(x(t,c0);A(t)), A(t0) = M

0 = F(c,x(t f ,c))

Khan and Barton (2014), Khan and Barton (2015), Pang and Stewart (2009).

20

Conclusions

◆ L-derivatives “computationally relevant”

generalized derivatives

◆ Can be computed automatically for broad

classes of functions

◆ Strong theory gives practically computable

generalized derivatives for:

Ø Implicit functions Ø Algorithms Ø ODE solutions Ø Linear programs Ø etc.

21

Acknowledgments

◆ Peter Stechlinski ◆ Novartis ◆ Statoil

22

Obstacles 3

◆ may be a strict subset of

∂f(x)

∂ fi(x)

i=1 m

∏

f :(x1, x2) ! x1+ | x2 | x1− | x2 | ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ ∂f(0) = 1 1 ⎡ ⎣ ⎢ 2s −1 1− 2s ⎤ ⎦ ⎥ : s ∈ 0,1 ⎡ ⎣ ⎤ ⎦ ⎧ ⎨ ⎩ ⎪ ⎫ ⎬ ⎭ ⎪ ∂ f1(0) × ∂ f2(0) = 1 1 ⎡ ⎣ ⎢ 2s1 −1 2s2 −1 ⎤ ⎦ ⎥ :(s1,s2) ∈ 0,1 ⎡ ⎣ ⎤ ⎦

2

⎧ ⎨ ⎪ ⎩ ⎪ ⎫ ⎬ ⎪ ⎭ ⎪

π 2∂f(0) = 2s −1 1− 2s ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ : s ∈ 0,1 ⎡ ⎣ ⎤ ⎦ ⎧ ⎨ ⎪ ⎩ ⎪ ⎫ ⎬ ⎪ ⎭ ⎪ π 2(∂f1(0) × ∂f2(0)) = 2s1 −1 2s2 −1 ⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ :(s1,s2) ∈ 0,1 ⎡ ⎣ ⎤ ⎦

2

⎧ ⎨ ⎪ ⎩ ⎪ ⎫ ⎬ ⎪ ⎭ ⎪

23

Lexicographic differentiation[1]

[1]: Y. Nesterov, Math. Program. B, 104 (2005) 669-700. ◆

is L-smooth at if it is locally Lipschitz continuous and directionally differentiable, and if, for any , the following functions exist:

◆

If the columns of span , then is linear

◆

Lexicographic subdifferential:

◆

The class of L-smooth functions is closed under composition, and includes all smooth functions and all convex functions f : X ∈ Rn → Rm x ∈ X M :=[m(1) ! m( p)] ∈ Rn×p fx,M

(0) :d ! f '(x;d)

fx,M

(1) :d ! [fx,M (0) ]'(m(1);d)

" fx,M

( p) :d ! [fx,M ( p−1)]'(m( p);d)

M Rn fx,M

( p)

∂Lf(x) ={JLf(x;M):M ∈Rn×n, det M ≠ 0} JLf(x;M):= Jf ( p)

x,M(0)

24

Inverse and implicit functions

1: Khan and Barton, submitted

◆ LD-derivatives for inverse functions:

Ø Suppose is L-smooth and locally invertible near , and is Lipschitz near Ø Result: is also L-smooth at Ø Result: is the unique solution of :

◆ LD-derivatives for implicit functions:

Ø Suppose is L-smooth, and there exists an implicit function such that Ø for each near Ø Result: is L-smooth at Ø Result: is the unique solution of:

f ˆ y, f(ˆ y) = ˆ z f −1 ˆ z f −1 ˆ z [f −1]'(ˆ z;M) N f '(ˆ y;N) = M h h(ˆ y, ˆ z) = 0, w h(w(z),z) = 0, z ˆ z ˆ z w w'(ˆ z;M) N h' ˆ y ˆ z ! " # # $ % & & ; N M ! " # $ % & ' ( ) ) * + , , = 0

25

Vector Forward AD Mode for LD-derivatives

◆

z

h(z) = f (g(z)) = f (z,z) = max{z,z,−z,−z}=| z |

f :!2 → ! :(x, y) " max{x, y,−x,− y} g :! → !2 : x " (x,x) h:! → ! : z " f #g Z f ={(x, y) ∈!2 : y = x or y = −x}

g(!) ⊂ Z f

26

Generalized Derivatives of Algorithms: MHEX model

F

i Ti in − Ti

• ut

( )

i∈H

∑

= f j ti

• ut − ti

in

( )

j∈ C

∑

min

p∈P EBP C p − EBP H p

( ) = 0

UA− ΔQk ΔTLM

k k∈K k≠|K|

∑

= 0 MHEX

F

1,T 1 in

F

1,T 1

• ut

f|C|,t|C|

• ut

f|C|,t|C|

• ut

! F

|H|,T|H|

• ut

F

|H|,T|H|

• ut

! ! ! f1,t1

in

f1,t1

• ut

27

◆ Consider the set of points which are either

kinks or endpoints on the composite curves

Ø Index this set of points with set Ø Each has an associated enthalpy value

T Q

k

Q

1 k

Q +

1 2

Q Q =

| | K

Q … …

K

k

Q k K ∈

Formulating a PC1 Area Constraint

28

◆ If both temperatures at adjacent points are

known, the interval can be treated as a two- stream heat exchanger

T Q

k

Q

1 k

Q +

1 2

Q Q =

| | K

Q … …

k

T

k

t

1 k

T

+ 1 k

t

+ k

Q Δ

ΔQk = UAkΔTLM

k

Formulating a PC1 Area Constraint

29

◆ Summing the area over all intervals gives the

total MHEX area: UA = ΔQk ΔTLM

k k∈K k≠|K|

∑

T Q

k

Q

1 k

Q +

1 2

Q Q =

k

Q Δ

| | K

Q … …

Formulating a PC1 Area Constraint

30

◆ Difficulties:

Ø The enthalpies and temperatures need to be sorted Ø Not all of the temperatures are known

◆ Consider (naïve) bubble sort: ◆ Only calculations involve taking min/max of two entries ◆ Same sequence and number of calculations regardless of

whether the input is well-sorted or not Naïve bubble sort is an composite PC1 function

Formulating a PC1 Area Constraint

31

◆ Finding the unknown temperatures involves solving one

• f these equations for :

◆ Easily solved using “if-else” logic

Ø Correctly calculates values of the unknown temperatures, but not their LD-derivatives ß not composite PC1

◆ Solution of the equation defines an implicit function

:

( ) ( )

• ut

in in

• ut

max{0, } max{0, } max{0, } max{0, } ( ( ) )

k k k i i i i H k k k j j j j C

Q F T T T T Q f t t t t

∈ ∈

− − − − = − − − − =

∑ ∑

y y

,

k k

T t

Calculating Unknown Temperatures

( )

, , ( )

k k

h T Q = y y η :!

ny × ! → !

( )

( ) , ,

k k

h T Q = y y

( )

( ) ( ) ,

k k

T Q η = y y y

( )

( )

ˆ ' ( ) ; ˆ '( ; ˆ ; ) ˆ ,

k k k

T Q Q η ⎡ ⎤ ⎛ ⎞ = ′⎜ ⎢ ⎜ ⎝ ⎣ ⎟ ⎠ ⎥ ⎦⎟ I y y I I y y

Generalization of the implicit function theorem

32

◆ Once all temperatures are found, the driving

force for each interval can be calculated

Ø Requires modification of standard log-mean temperature difference calculation: ◆ Finally, calculate the total area using:

( )

1 1 1

min min LM 1 1 1 1 2 1 ln( ) ln( )

max{ , } max{ , } if ,

• therwise

k k k k k k

k k k k k k k k T T k k k T T T T

T T T t T T T t T T T T T

+ + +

+ + + + Δ +Δ + Δ −Δ Δ − Δ

Δ = Δ − Δ = Δ − ⎧ Δ = Δ ⎪ Δ Δ Δ = ⎨ ⎪ ⎩

LM | | k k K k k K

Q UA T

∈ ≠

Δ = Δ

∑

Calculating Temperature Driving Force and Area

Continuously differentiable elemental function

33

◆ New extended pinch operator formulation: ◆ Relevant equation for the new formulation:

EBP

C p =

f j max 0,(T p − ΔTmin) − t j

in

{ }

(

j∈ C

∑

− max 0,(T p − ΔTmin) − t j

• ut

{ }

+max 0,(T p − ΔTmin) − tmax

{ }−max 0,tmin − (T p − ΔTmin) { })

EBP

H p =

Fi max 0,T p −Ti

• ut

{ }

(

i∈H

∑

− max 0,T p −Ti

in

{ }

−max 0,T min −T p

{ }+max 0,T p −T max { })

min

p∈P EBP C p − EBP H p

( ) = 0

Formulating the ΔTmin Equation

34

Complete Formulation

F

i Ti in − Ti

• ut

( )

i∈H

∑

− f j ti

• ut − ti

in

( )

j∈ C

∑

= 0 min

p∈P EBP C p − EBP H p

( ) = 0

UA− ΔQk ΔTLM

k k∈K k≠|K|

∑

= 0 Ti

• ut ≤ Ti

in,∀i ∈H

t j

• ut ≥ t j

in,∀j ∈C

3 equations (plus constraints on outlet temperatures) Solve for 3 unknown quantities using LP-Newton method (temperatures, flow rates, area, minimum approach temperature, etc.)

35

LNG process case study

◆ Two MHEX models plus three intermediate

compression/expansion operations

Ø 9 equations (4 nonsmooth) à solve for nine variables

36

LNG process data and variables

◆ Seven temperatures are taken as unknown in

each of the following cases

Ø 2 additional variables can be taken as unknowns (UA, approach temperature, more temperatures, etc.

37

LNG process example: Case I

min

4 K ? A T U Δ = =

min

4 K ? T UA = Δ =

◆ Given the minimum approach temperature,

calculate the area

38

LNG process example: Case I

HX-100: kW/K 97.54 UA = HX-101: kW/K 31.09 UA =

y1 y2 y3 y4 y5 y6 y7 365.07 K 225.44 K 193.35 K 95.08 K 180.85 K 357.14 K 95.14 K

39

Empirical Convergent Rate

20 40 60 80 100 120 140 Iteration

2 2 4 6 8 10

10 10 10 10 10 10 10

− − − − −

|| ( ) ||

k ∞

f y

40

LNG process example: Case II

min

? 120 (kW/K) UA T = = Δ

min

? 30 (kW/K) UA T = = Δ

◆ Given the area, calculate the minimum

approach temperature

41

LNG process example: Case II

min

HX-100: 2 K 2.6 T = Δ

min

HX-101: 6 K 1.2 T = Δ

y1 y2 y3 Y4 y5 y6 y7 365.07 K 217.66 K 196.09 K 95.08 K 174.75 K 362.46 K 91.85 K

42

LNG process example: Case III

min

4 K 85 (kW/K) UA T = Δ =

min

4 K 35 (kW/K) UA T = Δ =

◆ Given the area and the minimum approach

temperature, calculate (more) information about the streams

43

LNG process example: Case III

• ut

H2

199. K 06 T =

• ut

C3 = 168.

K 25 t

y1 y2 y3 y4 y5 y6 y7 365.07 K 231.10 K 195.06 K 95.08 K 194.58 K 352.12 K 97.52 K

44

PC RHS with non-PC Solution

◆

positive, independent of

! x = y, ! y = −x, ! z = | x |,

1 1

(0) (0) (0) x x y y z z = = = z(t) = z0 + | x0

t

∫

cos(s)+ y0 sin(s) |ds = z0 + |(rcos(θ))

t

∫

cos(s)+(rsin(θ))sin(s) |ds = z0 + x0

2 + y0 2

|cos(θ − s)

t

∫

|ds

zt(x0) for t ∈ [0,10] (y0 = 0, z0 = 0)

z2πk(x0, y0,0) = f (k) x0

2 + y0 2

⇒ z2π (x0, y0,0) x0 y0 θ

45

Existing dynamic sensitivities

◆

46

Existing dynamic sensitivities

◆

1 2 3 2 4 6 t Γ[x(2,⋅)](t)

1 2 3

• 5

5 t x(t,c)

47

Dynamic LD-derivatives

◆

48

Example 1

◆

0.5 1 1.5 2 2.5 3 1 2 3 4 5 6 Time (t) Generalized derivative bounds LNA bounds

• Lex. deriv. bounds

1 2 3

• 5

5 t x(t,c)

49

Singleton & Nonsingleton Trajectories

◆

OK OK Not OK