Brief Introduction to Numerical Optimization and its application to - PowerPoint PPT Presentation

Brief Introduction to Numerical Optimization and its application to Molecular Conformation Dongli Deng Bedˇ rich Soused´ ık (mentor) Supported in part by NSF award DMS-1521563 April 30, 2018

Introduction ◮ We studied four methods: ◮ Steepest Descend ◮ Conjugate gradient ◮ Newton’s method ◮ A quasi-Newton method ◮ We implemented the methods in MATLAB , and tested them on several small problems. ◮ We applied the methods to minimize Lennard-Jones potential.

Test function 1: (Ex. 2.1, Nocedal: Numerical Optimization, 2006) f ( x , y ) = 100( y − x 2 ) 2 + (1 − x ) 2 (minimum at (1 , 1)) − 400 xy + 400 x 3 − 2 + 2 x − 400 y + 1200 x 2 + 2 � � � � − 400 x ∇ f = , H = 200 y − 200 x 2 − 400 x 200 1 400 0.8 500 350 0.6 300 400 0.4 250 0.2 300 0 200 -0.2 200 150 -0.4 100 100 -0.6 -1 50 -0.8 0 0 1 0.5 -1 0 0 1 -0.5 -1 -1 -0.5 0 0.5 1

Test function 2: (Example 13.3, Sauer: Numerical analysis, Pearson, 2006) f ( x , y ) = 5 x 4 +4 x 2 − xy 3 +4 y 4 − x (minimum at (0 . 4923 , − 0 . 3643)) 20 x 3 + 8 xy − y 3 − 1 60 x 2 + 8 y − 3 y 2 � � � � ∇ f = , H = 4 x 2 − 3 y 2 x + 16 y 3 48 y 2 − 6 yx − 3 y 2 1 15 0.8 0.6 10 0.4 0.2 5 0 -0.2 -0.4 0 1 -0.6 0 -0.8 -5 -1 -0.5 -1 0 -1 0.5 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

Steepest Descend x k +1 = x k − sv s . . . step length (we can use Successive Parabolic Interpolation), v . . . direction of the steepest descend at x k (given by ∇ f ( x k ))

Algorithm of Steepest Descend 1: for n = 0 , 1 , 2 , . . . do v= ∇ f ( x i ) 2: Minimize f ( x − sv ) for scalar s = s ∗ 3: x i +1 = x i + s ∗ v 4: 5: end for

Steepest Descend for test problem 2 (Sauer) Step f ( x , y ) x y 1 1.0000 -1.0000 5.0000 5 0.1867 0.2136 -0.1444 10 0.3327 0.0418 -0.2530 15 0.4228 -0.2142 -0.4036 20 0.4877 -0.3561 -0.4573 25 0.4922 -0.3641 -0.4575 30 0.4923 -0.3643 -0.4575

Conjugate Gradient Method Consider minimizing, starting with quadratic function f ( x ) = 1 2 x T Ax − x T b ∇ f = Ax − b Finding the minimum is equivalent to solving Ax = b . One-dimensional line search: given search direction d , find step length α so that the function f ( x + α d ) is minimized 0 = ∇ f · d = ( α Ad − r ) T · d r T d r T r α = d T Ad = d T Ad r . . . residual of the linear system (given by −∇ f ( x ) = b − Ax ) .

0 = ∇ f · d = ( A ( x + α d ) − b ) · d 0 = ∇ f · d = ( Ax + α Ad − b ) · d Let r = b − Ax 0 = ( α Ad − r ) T · d α d T Ad − r T d = 0 r T d α = d T Ad

Algorithm of Conjugate Gradient Method 1: Let x 0 be the initial guess and set d 0 = r 0 = −∇ f . 2: for n = 0 , 1 , 2 , . . . do α i = α that minimizes f ( x i − 1 + α d i − 1 ) 3: x i = x i − 1 + α i d i − 1 4: r i = −∇ f ( x i ) 5: r T i r i β i = 6: r T i − 1 r i − 1 d i = r i + β i d i − 1 7: 8: end for

Conjugate Gradients for test problem 2 (Sauer) Step f ( x , y ) x y 1 0.0998 0.4114 0.0247 2 0.5372 -0.3240 -0.4324 3 0.5093 -0.3812 -0.4557 4 0.4944 -0.3776 -0.4569 5 0.4900 -0.3644 -0.4575 . . . 10 0.4923 -0.3643 -0.4575

Newton’s Method in one Dimension For solving f ( x ) = 0 we use x i +1 = x i − f ( x i ) f ′ ( x i ) , Applying the same idea to f ′ ( x ) = 0 gives x i +1 = x i − f ′ ( x i ) f ′′ ( x i ) . This can be equivalently written as f ′′ ( x i ) ( x i +1 − x i ) = − f ′ ( x i ) .

Newton’s Method in higher dimension Given an initial guess x 0 , for k = 0 , 1 , . . . , solve H ( x k ) v = −∇ f ( x k ) , update x k +1 = x k + v . Here � T � ∂ f ∂ f ∇ f = ∂ x 1 ( x 1 , ..., x n ) , ... ∂ x n ( x 1 , ..., x n ) , ∂ 2 f ( x 1 ,..., x n ) ∂ 2 f ( x 1 ,..., x n )   ... ∂ x 2 ∂ x 1 ∂ x n 1   ..     H f = .. .     ..     ∂ 2 f ( x 1 ,..., x n ) ∂ 2 f ( x 1 ,..., x n ) ... ∂ x 2 ∂ x 1 ∂ x n n

Quasi Newton’s method H . . . from now on, an approximation to the inverse of the Hessian , Update x as x k +1 = x k + α k p k α k . . . step length (from backtracking line search/Wolfe condition), = − H k ∇ f k (update of the search direction) , p k = x k +1 − x k (change of the displacement) , s k y k = ∇ f k +1 − ∇ f k (change of gradients of the functions) Update the (inverse of) the Hessian as H k +1 = ( I − ρ k s k y T k ) H k ( I − ρ k y k s T k ) + ρ k s k s T k , where 1 ρ k = k s k ) . ( y T

Algorithm of BFGS 1: Given starting point x 0 , convergence tolerance ǫ > 0 2: choose (inverse) Hessian approximation H 0 (commonly I ) 3: k ← 0 4: while | ∇ f k | > ǫ ; do p k = − H k ∇ f k 5: x k +1 = x k + α k p k 6: Compute H k +1 (updating H k ) 7: k ← k+1 8: 9: end while

Lennard-Jones Potential U ( r ) = 1 r 12 − 2 r . . . distance between two atoms r 6 , 1 0.5 0 -0.5 -1 0 0.5 1 1.5 2 2.5 3 r General form with n atoms n − 1 n � � 1 − 2 � � U ( x 1 , . . . , x n , y 1 , . . . , y n , z 1 , . . . , z n ) = , r 12 r 6 ij ij i =1 j = i +1 r ij . . . distance between atoms i and j

Gradient of the Lennard-Jones potential n − 1 n − 1 n n � � 1 − 2 � � � � U = U ij = r 12 r 6 ij ij i =1 j = i +1 i =1 j = i +1 and, for example, ∂ U ij = − 12( x i − x j ) + 12( x i − x j ) r 14 r 8 ∂ x i ij ij where � ( x i − x j ) 2 + ( y i − y j ) 2 + ( z i − z j ) 2 r ij =

Three atoms Left: initial guess (0 , − 5 , 0)(0 , 0 , 0)(0 , 5 , 0) (not global minimum) Right: initial guess (0 , 0 , 0)(0 , 0 , 2)(1 , 1 , 1) , comparison of fminunc and our implementation of BFGS 1 1.6 0.8 1.4 0.6 0.4 1.2 0.2 1 0 -0.2 0.8 -0.4 -0.6 0.6 -0.8 0.4 -1 1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 0 0.5 0.3 0.4 0.5 0.6 0.7 0.8 1 0.1 0.2

Four atoms Left: initial guess (0 , − 5 , 0)(0 , 0 , 0)(0 , 5 , 0)(0 , 10 , 0) (not global min.) Right: initial guess (0 , 0 , 0)(0 , 0 , 2)(1 , 1 , 1)(2 , 3 , 4), comparison of fminunc and our implementation of BFGS 1 2.2 0.8 2 0.6 1.8 0.4 1.6 1.4 0.2 1.2 0 1 -0.2 0.8 -0.4 0.6 -0.6 2 1 -0.8 0 -1 1.5 1 0.5 -1 4.5 4 3.5 3 2.5 2 1.5 1 0 -0.5

Six atoms 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 0.5 0.5 0 0 -0.5 -0.5 (We got the same result using fminunc in Matlab .)