numerical optimization and simulation 4311010
play

Numerical Optimization and Simulation (4311010) Manfred Gilli - PDF document

Dec 11, 2023 •206 likes •1.56k views

Numerical Optimization and Simulation (4311010) Manfred Gilli Department of Econometrics University of Geneva and Swiss Finance Institute Spring 2008 2 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  1. Iterative methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 General structure of algorithm for iterative methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 Solving least squares problems 61 Normal equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 Least squares solution with QR decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 Comparison between normal equations and QR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 Nonlinear equations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 Finding the roots of f ( x ) = 0 74 R´ esolution graphique. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 R´ esolution al´ eatoire. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 Recoupement par intervalles (bracketing) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 M´ ethode de la bisection. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 It´ erations de point fixe. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 M´ ethode de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 Systems of nonlinear equations 98 Iterative methods (Gauss-Seidel and Jacobi) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 Fix point method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 Newton method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 Quasi-Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 Damped Newton. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 Solution by minimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 Classical optimization paradigm 126 Definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 Classification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 More definitions and notations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 Unconstraint optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 Nonlinear least squares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 Direct search methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 Optimization heuristics 168 The standard optimization paradigm and its limits. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172 Heuristic paradigm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 189 Overview of optimization heuristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 190 Trajectory methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 194 Population based methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 199 Stochastics of solution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 207 Elements for classification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 211 Threshold accepting heuristic 219 Basic features. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 219 Local structure and neighborhood . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 221 Pseudo-code for TA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 224 Objective function. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 226 Constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 227 Neighborhood definition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 228 Threshold sequence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 229 TA with restarts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 231 The Shekel function 232 Searching the minimum of the Shekel function with TA 236 Neighborhood definition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 236 Computation of threshold sequence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 239 TA implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 240 2

  2. Visualization of working of TAH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 242 Genetic algorithm 244 Basic GA. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 244 An example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 245 Solving it with GA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 246 Generating the starting population . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 247 Selection of mating pool . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 248 Crossover. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 249 Mutation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 250 Selection of survivors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 251 Function GAH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 252 Differential Evolution algorithm (DE) 257 Introduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 257 Ackley’s function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 258 Initial population . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 259 Construction of new solutions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 261 Pseudo-code for DE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 266 Matlab code. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 267 Matlab code (cont’d) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 268 Matlab code (cont’d) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 269 Vectorized Matlab code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 270 Vectorized Matlab code (cont’d) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 271 Minimization of Ackley function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 272 3

  3. Outline Basics Computer arithmetic Representation of real numbers Machine precision Example of limitations of floating point arithmetic Numerical differenciation Forward-difference Backward-difference Central-difference Approximating second derivatives Partial derivatives How to choose h An example Measuring the error Absolute and relative error Combining absolute and relative error Numerical instability and ill conditioning Numerically unstable algorithm Numerically stable algorithm Ill conditioned problem Condition of a matrix Complexity of algorithms Order of complexity O ( · ) and classification Evolution of performance Solution of linear systems Ax = b Choice of method LU factorization Cholesky factorization QR facorization Sparse linear systems Sparse matrices in Matlab The slash operator in Matlab Practical considerations Iterative methods General structure of algorithm for iterative methods Solving least squares problems Normal equations Least squares solution with QR decomposition Comparison between normal equations and QR Nonlinear equations Finding the roots of f ( x ) = 0 R´ esolution graphique R´ esolution al´ eatoire Recoupement par intervalles (bracketing) M´ ethode de la bisection It´ erations de point fixe M´ ethode de Newton Systems of nonlinear equations Iterative methods (Gauss-Seidel and Jacobi) Fix point method Newton method Convergence 4

  4. Quasi-Newton Damped Newton Solution by minimization Classical optimization paradigm Definitions Classification More definitions and notations Unconstraint optimization Nonlinear least squares Direct search methods Optimization heuristics The standard optimization paradigm and its limits Heuristic paradigm Overview of optimization heuristics Trajectory methods Population based methods Stochastics of solution Elements for classification Threshold accepting heuristic Basic features Local structure and neighborhood Pseudo-code for TA Objective function Constraints Neighborhood definition Threshold sequence TA with restarts The Shekel function Searching the minimum of the Shekel function with TA Neighborhood definition Computation of threshold sequence TA implementation Visualization of working of TAH Genetic algorithm Basic GA An example Solving it with GA Generating the starting population Selection of mating pool Crossover Mutation Selection of survivors Function GAH Differential Evolution algorithm (DE) Introduction Ackley’s function Initial population Construction of new solutions Pseudo-code for DE Matlab code Matlab code (cont’d) Matlab code (cont’d) Vectorized Matlab code Vectorized Matlab code (cont’d) Minimization of Ackley function 5

  5. 2 Optimization: “The quest for the best” In practice model assumptions are rarely met we rather look for a “good” solution Numerical optimization: “All” the trouble comes from → In a computer it can happen that 1 + ǫ = 1 for ǫ � = 0 NOS (4311010) – M. Gilli Spring 2008 – 2 Introduction Optimization is a particularly broad and complex domain.This course aims at providing a structured overview of optimization problems and corresponding solution techniques. Emphasis will be on numerical and computational issues rather than on a formal presentation. We start with the basic standard techniques for unconstraint optimization of differentiable functions and recall related practical issues such as the solution of linear systems, the computation of derivatives etc. The presentation of some relevant optimization problems where classical methods fail to work will motivate the introduction of the heuristic optimization paradigm . NOS (4311010) – M. Gilli Spring 2008 – 3 Outline ◦ Basics. Computer arithmetic, numerical differentiation, solving linear systems, etc. ◦ Nonlinear equations. Equations in 1 dimension, systems of equations ◦ Overview of optimization problems and solution methods. Definitions, classification ◦ Unconstraint optimization. Gradient based methods, direct search ◦ Heuristic optimization paradigm. Overview and classification of optimization heuristics ◦ Local search methods. Simulated annealing, threshold accepting, tabu search ◦ Population based methods. Evolution based and genetic algorithms, ant systems and ant colony optimization, differential evolution, particle swarm optimization ◦ Monte Carlo methods. Random variable generation, quasi-Monte Carlo NOS (4311010) – M. Gilli Spring 2008 – 4 6

  6. Basics ◦ Approximations in numerical computation ◦ Computer arithmetic ◦ Rounding error ◦ Truncation error ◦ Numerical differentiation ◦ Finite difference approximation ◦ Automatic differentiation ◦ Solving linear systems ◦ LU decomposition ◦ Cholesky decomposition ◦ QR decomposition ◦ The slash operator in Matlab NOS (4311010) – M. Gilli Spring 2008 – 5 7

  7. Computer arithmetic 6 Two sources of errors appear systematically in numerical computation: ◦ Truncation errors due to the simplification of the mathematical model (e.g. replacement of a derivative by a finite difference, discretization); ◦ Rounding errors due to the fact that it is impossible to represent real numbers exactly in a computer. The support for all information in a computer is a “word” constituted by a sequence of bits (generally 32), a bit having value 0 or 1. 2 31 2 30 2 3 2 2 2 1 2 0 0 0 0 0 · · · · · · 0 0 1 0 1 0 This sequence of bits is then interpreted as the binary representation of an integer. As a consequence integers can be represented exactly in a computer. If all computations can be made with integers we face integer arithmetic . Such computations are exact. NOS (4311010) – M. Gilli Spring 2008 – 6 Representation of real numbers In a computer, a real variable x � = 0 is represented as x = ± n × b e n is the mantissa (or fractional part), b the base and e the exponent (base always 2). (Ex.: 12 . 153 = . 75956 × 2 4 or with base 10 we have 0 . 12153 × 10 2 ). In practice we partition the word in three parts, one containing e the second n and the first bit from left indicates the sign. · · · · · · · · · · · · ± � �� � � �� � e n Therefore, whatever size of the word, we dispose of a limited number of bits for the representation of the integers n et e . As a result it will be impossible to represent real numbers exactly. NOS (4311010) – M. Gilli Spring 2008 – 7 8

  8. Representation of real numbers To illustrate this fact, consider a word with 6 bits, with t = 3 the number of bits reserved for the mantissa, s = 2 bits reserved for the exponent and 1 bit for the sign. 0 0 0 1 0 0 4 0 1 1 1 0 1 5 e = n = ± 1 0 2 1 1 0 6 ���� � �� � e n 1 1 3 1 1 1 7 Normalizing the mantissa, i.e. imposing the first bit from left being 1, we have n ∈ { 4 , 5 , 6 , 7 } and defining the offset of the exponent as o = 2 s − 1 − 1 we have (0 − o ) ≤ e ≤ (2 s − 1 − o ) , i.e. e ∈ {− 1 , 0 , 1 , 2 } . The real number x then varies between (2 t − 1 ≤ n ≤ 2 t − 1) × 2 e NOS (4311010) – M. Gilli Spring 2008 – 8 Representation of real numbers Adding 1 to the right of (2 t − 1 ≤ n ≤ 2 t − 1) × 2 e we get n < 2 t and multiplying the whole expression by 2 − t we obtain the interval (2 − 1 ≤ n < 1) × 2 e − t for the representation of real numbers. NOS (4311010) – M. Gilli Spring 2008 – 9 Representation of real numbers Set of real numbers f = n × 2 e − t : 2 e − t n e = − 1 e = 0 e = 1 e = 2 1 / 16 1 / 8 1 / 4 1 / 2 1 / 4 1 / 2 1 2 =2 2 =4 1 0 0 =2 2 +2 0 =5 5 / 16 5 / 8 5 / 4 5 / 2 1 0 1 3 / 8 3 / 4 3 / 2 3 =2 2 +2 1 =6 1 1 0 7 / 16 7 / 8 7 / 4 7 / 2 =2 2 +2 1 +2 0 =7 1 1 1 We also observe that the representable real numbers are not equally spaced. 1 1 2 3 7 4 2 NOS (4311010) – M. Gilli Spring 2008 – 10 9

  9. Representation of real numbers Matlab uses double precision (64 bits) with t = 52 bits for the mantissa and e ∈ [ − 1023 , 1024] the range of integers for the exponent. We then have m ≤ f ≤ M with m ≈ 1 . 11 × 10 − 308 and M ≈ 1 . 67 × 10 308 . If the result of an expression is inferior to m we get an “ underflow ” if the result is superior to M we are in a situation of “ overflow ”. NOS (4311010) – M. Gilli Spring 2008 – 11 Representation of real numbers The notation float ( · ) represents the result of a floating point computation. A computer using the arithmetic of the preceding example ( t = 3 and s = 2 ) would produce with “ chopping ” the following numbers float (1 . 74) = 1 . 5 and float (0 . 74) = 0 . 625 float(0 . 74) float(1 . 74) 0 1 2 3 0 . 74 1 . 74 The result for float (1 . 74 − 0 . 74) = 0 . 875 , illustrates that even if the exact result is in F its representation might be different. NOS (4311010) – M. Gilli Spring 2008 – 12 Machine precision The precision of a floating point arithmetic is defined as the smallest positif number ǫ mach such that float (1 + ǫ mach ) > 1 The machine precision ǫ mach can be computed with the following algorithm 1: e = 1 2: while 1 + e > 1 do 3: e = e/ 2 4: end while 5: ǫ mach = 2 e With Matlab we have ǫ mach ≈ 2 . 2 × 10 − 16 . NOS (4311010) – M. Gilli Spring 2008 – 13 Rounding errors and precision Floating point computation is not associative. We can observe that for e = ǫ mach / 2 we obtain the following results: float ((1 + e ) + e ) = 1 float (1 + ( e + e )) > 1 NOS (4311010) – M. Gilli Spring 2008 – 14 10

  10. Example of limitations of floating point arithmetic Sometimes the approximation obtained with floating point arithmetic can be surprising. Consider the expression ∞ 1 � n = ∞ n =1 Computing this sum will, of course, produce a finite result. What might be surprising is that this sum is certainly inferior to 100. The problem is not generated by an “ underflow ” of 1 /n or an “ overflow ” of the partial sum � k n =1 1 /n but by the fact that for n verifying 1 /n < ǫ mach n − 1 1 � k k =1 the sum remains constant. NOS (4311010) – M. Gilli Spring 2008 – 15 Example of limitations of floating point arithmetic (contd.) To show this we consider e < ǫ mach and we can write float ( 1 + e ) = 1 � n − 1 � n − 1 1 1 1 float ( + ) = k =1 k =1 k n k 1 /n float ( 1 + ) = 1 � n − 1 1 k =1 k We can easily experiment that � n − 1 1 k < 100 for values of n which can be considered in practice k =1 100 50 0 0 2 4 6 8 10 6 x 10 and given that ǫ mach ≈ 2 × 10 − 16 we establish 1 /n 100 = 2 × 10 − 16 from which we conclude that n is of order 2 × 10 14 . For a computer with a performance of 100 Mflops, the sum will stop to increase after (2 × 2 × 10 14 ) / 10 8 ) = 4 × 10 6 seconds, i.e. 46 days of computing without exceeding the value of 100. NOS (4311010) – M. Gilli Spring 2008 – 16 11

  11. Numerical differenciation 17 Approximation of derivatives We want to replace derivatives (partial) by numerical approximations. Given f : R → R a continuous function with f ′ continue. For x ∈ R we can write: f ( x + h ) − f ( x ) f ′ ( x ) = lim (1) h h → 0 If, instead of letting h go to zero, we consider “small” values for h , we get an approximation for f ′ ( x ) Question: How to choose h to get a good approximation if we work with floating point arithmetic? NOS (4311010) – M. Gilli Spring 2008 – 17 Approximation of derivatives (cont.d) We consider f with derivatives of order n + 1 For h finite, Taylor expansion writes f ( x + h ) = f ( x ) + h f ′ ( x ) + h 2 2 f ′′ ( x ) + h 3 6 f ′′′ ( x ) + · · · + h n n ! f ( n ) ( x ) + R n ( x + h ) with the remainder h n +1 ( n + 1)! f ( n +1) ( ξ ) R n ( x + h ) = avec ξ ∈ [ x, x + h ] ξ not known !! Thus Taylor series allow to approximate a function with a degree of precision which equals the remainder. NOS (4311010) – M. Gilli Spring 2008 – 18 Forward-difference Consider Taylors expansion up to order two f ( x + h ) = f ( x ) + h f ′ ( x ) + h 2 2 f ′′ ( ξ ) avec ξ ∈ [ x, x + h ] from which we get f ′ ( x ) = f ( x + h ) − f ( x ) − h 2 f ′′ ( ξ ) (2) h Expression (2) decomposes f ′ ( x ) in two parts, the approximation of the derivative and the truncation error This approximation is called forward-difference . NOS (4311010) – M. Gilli Spring 2008 – 19 12

  12. Backward-difference Choosing h negativewe get f ′ ( x ) = f ( x ) − f ( x − h ) + h 2 f ′′ ( ξ ) (3) h This defines the backward-difference approximation. NOS (4311010) – M. Gilli Spring 2008 – 20 Central-difference Consider the difference f ( x + h ) − f ( x − h ) : f ( x + h ) = f ( x )+ h f ′ ( x ) + h 2 2 f ′′ ( x )+ h 3 6 f ′′′ ( ξ + ) avec ξ + ∈ [ x, x + h ] f ( x − h ) = f ( x ) − h f ′ ( x ) + h 2 2 f ′′ ( x ) − h 3 6 f ′′′ ( ξ − ) avec ξ − ∈ [ x − h, x ] We have f ( x + h ) − f ( x − h ) = 2 h f ′ ( x ) + h 3 � � f ′′′ ( ξ + ) + f ′′′ ( ξ − ) 6 If f ′′′ is continuous we can consider the mean f ′′′ ( ξ ) , ξ ∈ [ x − h, x + h ] h 3 = 2 h 3 f ′′′ ( ξ + ) + f ′′′ ( ξ − ) � � f ′′′ ( ξ + ) + f ′′′ ( ξ − ) 3 3 2 � �� � f ′′′ ( ξ ) − h 2 f ′ ( x ) = f ( x + h ) − f ( x − h ) 3 f ′′′ ( ξ ) 2 h Defines the approximation by central-difference More precise as truncation error is of order O ( h 2 ) NOS (4311010) – M. Gilli Spring 2008 – 21 Approximating second derivatives Developing the sum f ( x + h ) + f ( x − h ) f ( x + h ) = f ( x )+ h f ′ ( x ) + h 2 2 f ′′ ( x )+ h 3 6 f ′′′ ( ξ + ) avec ξ + ∈ [ x, x + h ] f ( x − h ) = f ( x ) − h f ′ ( x ) + h 2 2 f ′′ ( x ) − h 3 6 f ′′′ ( ξ − ) avec ξ − ∈ [ x − h, x ] we get − h 2 f ′′ ( x ) = f ( x + h ) − 2 f ( x ) + f ( x − h ) 3 f ′′′ ( ξ ) avec ξ ∈ [ x − h, x + h ] h 2 There exist several different approximation schemes, also of higher order (see (Dennis and Schnabel, 1983)) NOS (4311010) – M. Gilli Spring 2008 – 22 13

  13. Partial derivatives Given a function of two variables f ( x, y ) then the approximation, using central-difference, of the partial derivative with respect to x is f x = f ( x + h x , y ) − f ( x − h x , y ) 2 h x Approximating with a central-difference the partial derivative of f x with respect to y gives f ( x + h x ,y + h y ) − f ( x − h x ,y + h y ) − f ( x + h x ,y − h y ) − f ( x − h x ,y − h y ) 2 h x 2 h x f xy = 2 h y 1 = ( f ( x + h x , y + h y ) − f ( x − h x , y + h y ) 4 h x h y − f ( x + h x , y − h y ) + f ( x − h x , y − h y )) NOS (4311010) – M. Gilli Spring 2008 – 23 How to choose h How to choose h if we use floating point arithmetic ? To reduce the truncation error we would be tempted to choose h as small as possible. However we also must consider the rounding error !! If h is too small, we have float ( x + h ) = float ( x ) and even if float ( x + h ) � = float ( x ) , � � � � we can have float f ( x + h ) = float f ( x ) if the function varies very slowly NOS (4311010) – M. Gilli Spring 2008 – 24 Truncation error for forward-difference − h 2 f ′′ ( ξ ) Denote h ( x ) = f ( x + h ) − f ( x ) f ′ h the approximation of the derivative corresponding to a given value of h If we accept that in the truncation error − h 2 f ′′ ( ξ ) we have | f ′′ ( t ) | ≤ M for all t in the domain of our application we get h ( x ) − f ′ ( x ) | ≤ h | f ′ 2 M bound for the truncation error, indicates that we can reach any precision given we choose h sufficiently small. However, in floating point arithmetic it is not possible to evaluate f ′ h ( x ) exactly (rounding errors !!) NOS (4311010) – M. Gilli Spring 2008 – 25 14

  14. How to choose h (cont’d) Consider that the rounding error for evaluating f ( x ) has the following bound � � � � � � f ( x ) − f ( x ) � ≤ ǫ � float then the rounding error (for finite precision computations) for an approximation of type forward-difference is bounded by � � � � � ≤ 2 ǫ f ′ − f ′ � � � float h ( x ) h ( x ) h � � � f ( x + h ) − f ( x ) � f ′ h ( x ) = float float h � � � � − float float f ( x + h ) f ( x ) = h ǫ + ǫ h Finally, the upper bound for the sum of both sources of errors is � � � � � ≤ 2 ǫ + M � f ′ − f ′ ( x ) � h ( x ) 2 h � float h Necessitates compromise between reduction of the truncation error and reduction of rounding error Consequence: The precision we can achieve is limited NOS (4311010) – M. Gilli Spring 2008 – 26 h is chosen as power of 2 and therefore can be represented without rounding error NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 26 15

  15. How to choose h (cont’d) Minimization problem of g ( h ) = 2 ǫ h + h 2 M h g ( h ) ⇒ g ′ ( h ) = 0 min − 2 ǫ h 2 + M = 0 2 4 ǫ h 2 = M � ǫ h = 2 M Bound reaches its minimum for � ǫ h = 2 (4) M In practice ǫ corresponds to the machine precision ǫ mach and if we admit that M is of order one, thus we set M = 1 and h = 2 √ ǫ mach With Matlab we have ǫ mach ≈ 2 × 10 − 16 and h ≈ 10 − 8 NOS (4311010) – M. Gilli Spring 2008 – 27 How to choose h (cont’d) h + h 2 Considering central-differences we have g ( h ) = 2 ǫ 3 M h g ( h ) ⇒ g ′ ( h ) = 0 min − 2 ǫ h 2 + 2 Mh = 0 3 6 ǫ h 3 = 2 M � 3 ǫ � 1 / 3 h = M h ≈ 10 − 5 NOS (4311010) – M. Gilli Spring 2008 – 28 16

  16. An example To illustrate the influence of the choice of h on the precision of the approximation of the derivative, consider the function f ( x ) = cos( x x ) − sin( e x ) and its derivative f ′ ( x ) = − sin( x x ) x x � � − cos( e x ) e x . log( x ) + 1 NOS (4311010) – M. Gilli Spring 2008 – 29 How to choose h (cont’d) f(x) = cos(x x ) − sin(e x ) 0 1.5 10 1 0.5 −5 10 Erreur rélative 0 −0.5 −10 10 −1 −1.5 −15 −2 10 −15 −10 −5 0 0.5 1 1.5 2 2.5 10 10 10 10 x h Relative error for the approximation of the derivative for x = 1 . 77 , as a function of h in the range of 2 0 et 2 − 52 , corresponding to ǫ mach for Matlab. NOS (4311010) – M. Gilli Spring 2008 – 30 17

  17. Measuring the error 31 Absolute and relative error Measuring the error e between an approximated value ˆ x and an exact value x . ◦ Absolute error | ˆ x − x | For ˆ x = 3 and x = 2 absolute error is 1, cannot be considered as “small” x = 10 9 + 1 and x = 10 9 can be considered “small” with respect to x For ˆ ◦ Relative error | ˆ x − x | Defined if x � = 0 not defined if x = 0 | x | If x = 0 the absolute error would do the job !! The relative error for the previous example is 0 . 5 and 10 − 9 respectively (“small” for the second case) NOS (4311010) – M. Gilli Spring 2008 – 31 Combining absolute and relative error Combination between absolute and relative error | ˆ x − x | | x | + 1 ◦ Avoids problems if x is near zero ◦ Has properties of the relative error if | x | ≫ 1 ◦ Has properties of absolute error if | x | ≪ 1 NOS (4311010) – M. Gilli Spring 2008 – 32 18

  18. Numerical instability and ill conditioning 33 Numerical instability If the “quality” of a solution is not acceptable it is important to distinguish between the following two situations: ◦ Errors are amplified by the algorithm → numerical instability ◦ Small perturbations of data generate large changes in the solution → ill conditioned problem NOS (4311010) – M. Gilli Spring 2008 – 33 Numerically unstable algorithm Example of a numerically unstable algorithm: solution of ax 2 + bx + c = 0 Analytical solutions are: √ √ b 2 − 4 ac b 2 − 4 ac x 1 = − b − x 2 = − b + . 2 a 2 a Algorithm: Transcription of analytical solution √ b 2 − 4 ac 1: ∆ = 2: x 1 = ( − b − ∆) / (2 a ) 3: x 2 = ( − b + ∆) / (2 a ) For a = 1 , c = 2 and a floating point arithmetic with 5 digits of precision the algorithm produces the following solutions for different values of b : | float( x 2 ) − x 2 | b ∆ float(∆) float( x 2 ) x 2 | x 2 | +1 1 . 55 × 10 − 6 5 . 2123 4 . 378135 4 . 3781 − 0 . 41708 − 0 . 4170822 1 . 47 × 10 − 3 121 . 23 121 . 197 121 . 20 − 0 . 01500 − 0 . 0164998 1212 . 3 1212 . 2967 1212 . 3 0 − 0 . 001649757 Catastrophic cancelation Attention: float ( x 2 ) = ( − b + float ( ∆)) / (2 a ) NOS (4311010) – M. Gilli Spring 2008 – 34 Numerically stable algorithm Alternative algorithm, exploiting x 1 x 2 = c/a (Vi` ete theorem): √ b 2 − 4 ac 1: ∆ = 2: if b < 0 then x 1 = ( − b + ∆) / (2 a ) 3: 4: else x 1 = ( − b − ∆) / (2 a ) 5: 6: end if 7: x 2 = c/ ( a x 1 ) b float(∆) float( x 1 ) float( x 2 ) x 2 1212 . 3 1212 . 3 − 1212 . 3 − 0 . 0016498 − 0 . 001649757 Avoids catastrophic cancelation (loss of significant digits when subtracting two large numbers). NOS (4311010) – M. Gilli Spring 2008 – 35 19

  19. Ill conditioned problem Measures the sensitivity of the solution of a linear system Ax = b with respect to a perturbation of the elements of A Example: � . 780 � . 217 � � . 563 A = b = . 913 . 659 . 254 The exact solution is x = [1 − 1] ′ (generally not known !!). Matlab computes � � . 99999999991008 x = A \ b = − . 99999999987542 This is an ill conditioned problem !!! NOS (4311010) – M. Gilli Spring 2008 – 36 20

  20. Ill conditioned problem Consider the matrix � � . 001 . 001 E = − . 002 − . 001 and the perturbed system ( A + E ) x E = b The new solution is � � − 5 . 0000 x E = 7 . 3085 5 5 0 0 −5 −5 −5 0 5 −5 0 5 5 4.5419 4.5419 4.5419 4.5418 4.5417 −1 4.5417 4.5416 4.5416 −5 −5 1 5 −3.0001−3.0001 −3 −3 −2.9999−2.9999 NOS (4311010) – M. Gilli Spring 2008 – 37 Condition of a matrix The condition of a matrix κ ( A ) is defined as κ ( A ) = � A − 1 � � A � Matlab function cond computes the condition of a matrix. For our example we have cond ( A ) = 2 . 2 × 10 6 If cond ( A ) > 1 / sqrt ( eps ) , we should worry about our numerical results !! For eps = 2 . 2 × 10 − 16 , we have 1 / sqrt ( eps ) = 6 . 7 × 10 8 (We loose half of the significant digits !!!) NOS (4311010) – M. Gilli Spring 2008 – 38 21

  21. Complexity of algorithms 39 ◦ Many different algorithms to solve same problem ◦ How to compare them in order to choose the best Precise comparison generally very difficult → use notion of complexity ◦ Problem size → number of elements to be processed e.g. linear algebra algorithms: Ax = b order n of A if A ∈ R n × m size given by n and m graph theory: G = ( V, E ) , n = | V | and m = | E | Time necessary to execute the algorithm expressed as a function of problem size → independent from executing platform ◦ Operation count (flop): only 4 elementary operations + − × \ considered ◦ Only order of function counting operations considered 3 + n 2 + 2 e.g. Complexity of algorithm with n 3 3 n elementary operations is of order O ( n 3 ) NOS (4311010) – M. Gilli Spring 2008 – 39 Order of complexity O ( · ) and classification Definition: Function g ( n ) is O ( f ( n )) if there exist constants c o and n o such that g ( n ) is inferior to c o f ( n ) for all n > n 0 g ( n ) = (5 / 4) n 3 + (1 / 5) n 2 + 500 O ( n 3 ) g ( n ) is as c 0 f ( n ) > g ( n ) for n > n 0 c 0 = 2 and n 0 = 9 g(n) c 0 f(n) f(n)=n 3 1374 6 8.8256 12 15 x 10 11 g(n) f(n)=n 3 10 n 2 5 0 0 2000 4000 6000 8000 10000 Algorithms are classified into 2 groups (according to the order of the function counting elementary operations): ◦ polynomial ◦ non-polynomial Only algorithms from the polynomial class can be considered as efficient !! Performance of computer measured in Mflops ( 10 6 flops per second) NOS (4311010) – M. Gilli Spring 2008 – 40 22

  22. Evolution of performance 4 10 2 10 Machine Mflops Year Mflops 8086/87 ` a 8 MHz 0.03 1985 Cray X–MP 33 1985 0 Pentium 133 MHz 12 1995 10 Pentium II 233 MHz 28 1997 Pentium III 500 MHz 81 1999 Pentium IV 2.8 GHz 1000 2003 Intel U2500 1.2 GHz 1000 2006 −2 10 85 95 97 99 03 06 Earth Simulator 30 TFlops 2003 www.top500.org Fantastic increase of cheap computing power !!! We will make use of it !!! NOS (4311010) – M. Gilli Spring 2008 – 41 23

  23. Solution of linear systems Ax = b 42 Solution of linear systems Ax = b Find solution to the following system of equations       a 11 x 1 + a 12 x 2 + · · · + a 1 n x n = b 1 a 11 a 12 · · · a 1 n x 1 b 1 a 21 x 1 + a 22 x 2 + · · · + a 2 n x n = b 2 a 21 a 22 · · · a 2 n x 2 b 2             ≡ = . . . . . . . . . ... . . . .  . . .   .   .  . . . . . . . . .       a n 1 x 1 + a n 2 x 2 + · · · + a nn x n = b n a n 1 a n 2 · · · a nn x n b n � �� � � �� � � �� � x A b ◦ Solution exists and is unique if and only if | A | � = 0 The notation for the solution is x = A − 1 b , might suggest that we need to compute A − 1 ◦ An efficient numerical method never computes A − 1 to solve the system nor | A | to check whether a ◦ solution exists !!! NOS (4311010) – M. Gilli Spring 2008 – 42 Choice of method Choice of appropriate method depends on structure and particular quantification of matrix A • Properties of structure ◦ sparse ◦ concentration of nonzero elements near diagonal ( a ij = 0 if | i − j | > b , b is called the bandwidth) diagonal ( b = 0 ), tridiagonal matrix ( b = 1 ) ◦ triangular matrix, block structure • Properties of quantification: symmetric ( A = A ′ ) ◦ positive definite ( x ′ Ax > 0 , ◦ ∀ x � = 0 ) Two broad categories of methods: ◦ Direct methods (resort to factorization) ◦ Iterative methods (stationary, non-stationary) NOS (4311010) – M. Gilli Spring 2008 – 43 LU factorization Direct methods resort to a factorization of matrix A (transform the system into one easier to solve) LU factorization (method of choice if A has no particular structure nor quantification) U = A L L Lower triangular matrix U Upper triangular matrix 3 n 3 − 1 2 n 2 − 1 2 Operation count: 6 n + 1 NOS (4311010) – M. Gilli Spring 2008 – 44 24

  24. I llustrate complexity with example n 3 + n 2 + n for n = 100 , 1 000 , 10 000 NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 44 Solution of Ax = b In Ax = b we replace A by its factorization U x = b L y y = Solve Ly = b ( forward substitution ) b L U y x = Solve Ux = y ( backward substitution ) NOS (4311010) – M. Gilli Spring 2008 – 45 LU factorization with Matlab ◦ Matlab uses the syntax x = A \ b to solve Ax = b (no matter what A is) ◦ LU factorization Matlab: [L,U] = lu(A); But L is a permuted triangular matrix !!! We can not apply a forward substitution ◦ [L,U,P] = lu(A); produces a triangular matrix L P is a permutation matrix such that LU = PA [L,U,P] = lu(A); x = U \ ( L \ P*b ); If the system is solved only once, we simply code x = A \ b (see ExSL0.m ) NOS (4311010) – M. Gilli Spring 2008 – 46 25

  25. ExSL0.m Generate A with Matlab and compute LU n = 5; x = ones(n ,1); A = unidrnd (40,n,n) b = A*x; [L,U] = lu(A) x = L\(U\b) [L,U,P] = lu(A) x = U\(L\P*b) NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 46 Cholesky factorization If A is symmetric and definite positive ( A = X ′ X with X full column rank) we use the Cholesky factorization R = A R ′ 3 n 3 + n 2 + 8 1 ◦ Operation count: 3 n − 2 flops ◦ Matlab: [R,p] = chol(A); ◦ Cholesky decomposition used to test numerically whether A is definite positive. If p � = 0 only the submatrix A(1:p,1:p) is definite positive ◦ R is used to generate multivariate random variables x with given Σ E ( xx ′ ) = E ( R ′ uu ′ x = R ′ u R ) = R ′ R = Σ u ∼ N(0 , I ) and ���� I NOS (4311010) – M. Gilli Spring 2008 – 47 Example: ExSL0.m Q = [1.00 0.30 -.30 0.30 2.00 0.70 -.30 0.70 1.50]; [R,p] = chol(Q); if p ~= 0, error(’Q not definite positive ’); end % u = randn (3 ,1); z = R’*u % n = 5000; U = randn (3,n); Z = R’*U; cov(Z’) NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 47 26

  26. QR facorization If A is square or rectangular and has full column rank ( A ∈ R m × n ) we compute the QR decomposition R = Q 1 Q 2 A Q is orthogonal, i.e. Q ′ Q = I ◦ Complexity: 4 m 2 n ◦ ◦ Matlab: [Q,R] = qr(A); NOS (4311010) – M. Gilli Spring 2008 – 48 Solution of Ax = b Replace A by its factorization R Q x = b R Q ′ y Solve Rx = Q ′ b x = ( backward substitution ) [Q,R] = qr(A); x = R \ (Q’*b); NOS (4311010) – M. Gilli Spring 2008 – 49 Sparse linear systems • Sparse matrix: any matrix with a sufficient number of zero elements such that there is a gain in considering them (avoid redundant operations and storage) • In many practical problems in economics and finance arise very large but sparse linear systems (e.g. discretization of PDEs). Cannot be solved without resorting to appropriate algorithms. • Matlab has sparse direct methods NOS (4311010) – M. Gilli Spring 2008 – 50 27

  27. Tridiagonal systems       1 u 1 r 1 d 1 q 1 l 1 1 u 2 r 2 p 1 d 2 q 2             l 2 1 u 3 r 3 = p 2 d 3 q 3             l 3 1 u 4 r 4 p 3 d 4 q 4       l 4 1 u 5 p 4 d 5 � �� � � �� � � �� � L U A Factorization of tridiagonal A with p i , i = 1 , . . . , n − 1 the lower diagonal, d i , i = 1 , . . . , n the diagonal and q i , i = 1 , . . . , n − 1 the upper diagonal. We verify that r i = q i , i = 1 , . . . , n − 1 . 1: u 1 = d 1 2: for i = 2 : n do 3: l i − 1 = p i − 1 /u i − 1 4: u i = d i − l i − 1 q i − 1 5: end for Operation count: 4 n − 4 flops NOS (4311010) – M. Gilli Spring 2008 – 51 Irregular sparse matrices Different storage schemes exist (problem dependent, Matlab → column wise) For A with n = 3 and nnz = 4 (nonzeros) in Matlab we have: Columns   1 2 3 0 1 0 4 × n (+1) Column pointers: 1 3 A = 2 0 3   4 0 0 1 2 3 4 2 3 1 2 nnz Row index: nnz 2 4 1 3 Elements: This storage scheme needs n integers for the row pointers, nnz integers for the row indices and nnz reals for the elements. An integer needs 4 bytes and a real 8 bytes, total bytes for the storage of a matrix is 12 nnz + 4 n bytes NOS (4311010) – M. Gilli Spring 2008 – 52 28

  28. Sparse matrices in Matlab • Conversion to sparse matrix not automatic • User must execute an explicit command • Once initialized, the sparse storage is propagated (an operation with a sparse matrix produces a sparse result, except for addition and substraction) B = sparse(A); C = full(B); • Practically sparse matrices are created directly: S = sparse(i,j,s,m,n,nzmax); Creates a sparse matrix of size m × n with S i ( k ) ,j ( k ) = s ( k ) S = sparse(i,j,s,m,n); % nzmax = length(s) S = sparse(i,j,s); % m = max(i) n = max(j) S = sparse(m,n); % S = sparse([] ,[] ,[] ,n,m) S = sparse([2 3 1 2] ,[1 1 2 3] ,[2 4 1 3]) NOS (4311010) – M. Gilli Spring 2008 – 53 Functions for sparse matrices in Matlab help sparfun speye(n,m) sprandn (n,m,d) spdiags spalloc (n,m,nzmax) issparse(S) spfun % C = spfun(’exp ’,B) spy(S) condest (S) eigs sprank [p,q,r,s] = dmperm(S) Ex.: /Teaching/MNE/Ex/ExSparse.m NOS (4311010) – M. Gilli Spring 2008 – 54 29

  29. The slash operator in Matlab Pseudo-code for how \ works with full matrices: if size(A,1) == size(A ,2) % A is square if isequal (A,tril(A)) % A is lower triangular x = A \ b; % Forward substitution on b elseif isequal (A,triu(A)) % A is upper triangular x = A \ b; % Backward substitution on b else if isequal (A,A’) % A is symmetric [R,p] = chol(A); if (p == 0) % A is symmetric positive definite x = R \ (R’ \ b); % Forward and backward substitution return end end [L,U,P] = lu(A); % general , square A x = U \ (L \ (P*b)); % Forward and backward substitution end else % A is rectangular [Q,R] = qr(A); x = R \ (Q’ * b); end NOS (4311010) – M. Gilli Spring 2008 – 55 Practical considerations • If same linear system Ax = b has to be solved for different b (stored in n columns of B , i.e. Ax = B ) we proceed as: [L,U,P] = lu(A); for j = 1:n x = U \ ( L \ P*B(:,j) ); end Matrix is factorized only once and in the following we only execute forward and backward substitutions (we simply code X = A \ B ) impossible to beat the backslash operator !!! If B = I , X = A − 1 How to compute s = q ′ A − 1 r without forming A − 1 ? • A − 1 r is the solution of Ax = r , → s = q’ * (A\r) NOS (4311010) – M. Gilli Spring 2008 – 56 Iterative methods • Stationary (Jacobi, Gauss-Seidel, SOR) • Non stationary (Krylov methods, QMRES, MINRES, BiCGstab, etc.) • Direct methods produce an exact solution in a number of finite operations Indirect methods produce an sequence { x ( k ) } , k = 0 , 1 , 2 , . . . • Produce an approximation of the solution in a number of infinite operations (in practice approximation acceptable in a limited number of operations) • Iterative methods are easy to implement (only matrix vector products) Makes it easy to exploit the sparse structure of a problem • However, efficiency depends on the speed of convergence !!! ( generally not guaranteed ) NOS (4311010) – M. Gilli Spring 2008 – 57 30

  30. Stationary iterative methods a 11 x 1 + a 12 x 2 + a 13 x 3 = b 1 a 21 x 1 + a 22 x 2 + a 23 x 3 = b 2 a 31 x 1 + a 32 x 2 + a 33 x 3 = b 3 We isolate the diagonal element in each row: x 1 = ( b 1 − a 12 x 2 − a 13 x 3 ) /a 11 x 2 = ( b 2 − a 21 x 1 − a 23 x 3 ) /a 22 x 3 = ( b 3 − a 31 x 1 − a 32 x 2 ) /a 33 x i ’s on the right-hand side unknown! Consider x ( k ) an approximation for x = A − 1 b , then we can compute a new approximation as x ( k +1) = ( b 1 − a 12 x ( k ) − a 13 x ( k ) 3 ) /a 11 1 2 x ( k +1) = ( b 2 − a 21 x ( k ) − a 23 x ( k ) 3 ) /a 22 2 1 x ( k +1) = ( b 3 − a 31 x ( k ) − a 32 x ( k ) 2 ) /a 33 3 1 This defines the Jacobi iteration NOS (4311010) – M. Gilli Spring 2008 – 58 Stationary iterative methods Use new approximations in the right-hand side as soon as they are available x ( k +1) =( b 1 − a 12 x ( k ) − a 13 x ( k ) 3 ) /a 11 1 2 x ( k +1) =( b 2 − a 21 x ( k +1) − a 23 x ( k ) 3 ) /a 22 2 1 x ( k +1) =( b 3 − a 31 x ( k ) − a 32 x ( k +1) ) /a 33 3 1 2 This defines the Gauss-Seidel iteration In practice we overwrite the x vector: Give initial solution x (0) ∈ R n i − 1 i i + 1 1 2 n for k = 0 , 1 , . . . until convergence do x : for i = 1 : n do � � b i − � x i = j � = i a ij x j /a ii x ( k +1) x ( k ) end for x ( k +1) end for i NOS (4311010) – M. Gilli Spring 2008 – 59 31

  31. General structure of algorithm for iterative methods Initialize x (1) , x (0) , ε and maxit while ˜ converged( x (0) , x (1) , ε ) do x (0) = x (1) (Store previous iteration in x (0) ) Compute x (1) with Jacobi, Gauss-Seidel or SOR Stop if number of iterations exceeds maxit end while function res = converged (x0 ,x1 ,tol) res = all(abs(x1 -x0 )./( abs(x0 )+1) < tol ); ◦ Convergence of Jacobi is sensitive to the normalization of the equations ◦ Convergence of Gauss-Seidel is sensitive to the normalization and the ordering of the equations No practical method for checking convergence properties exist! NOS (4311010) – M. Gilli Spring 2008 – 60 32

  32. Solving least squares problems 61 We consider the overdetermined system Ax ≈ b where b are the observations, A the independent variables and x the parameters. Residuals are r = b − Ax (Corresponding notation in econometrics: y ≡ b , X ≡ A and β ≡ x ) Least squares minimizes � b − Ax � 2 ◦ 2 (other objective functions can be considered !!!) This minimization is computationally very efficient ( reason for choice !! ) 2 at the solution in order to compute σ 2 = � b − Ax � 2 Evaluate � b − Ax � 2 ◦ 2 / ( m − n ) Compute efficiently ( A ′ A ) − 1 (or a subset of elements) ◦ Available methods: ◦ Normal equations ◦ QR decomposition ◦ SVD NOS (4311010) – M. Gilli Spring 2008 – 61 Normal equations 1 2 r ( x ) ′ r ( x ) min x ∈ R n The first order conditions are 2 A ′ Ax − 2 A ′ b = 0 from where we get the normal equations A ′ Ax = A ′ b Can also be obtained from the orthogonality condition A ′ r = A ′ ( b − Ax ) = 0 A ′ Ax = A ′ b → If A has full column rank A ′ A is positive definite and the linear system forming the normal equations can be solved using Cholesky factorization NOS (4311010) – M. Gilli Spring 2008 – 62 33

  33. Normal equations 1: Compute C = A ′ A � � A ′ b C 2: Form C = b ′ A b ′ b � G � 0 ′ 3: Compute G = ( C = G G Cholesky decomposition) z ′ ρ 4: Solve G ′ x = z (Triangular system) 5: σ 2 = ρ 2 / ( m − n ) 6: Compute T = G − 1 (Triangular matrix) 7: S = σ 2 T ′ T (Variance-covariance matrix) S = ( A ′ A ) − 1 � 1 /g ii i = j T = G − 1 = ( GG ′ ) − 1 t ij = � � i − 1 � − k = j g ik t kj /g ii i � = j = G ′− 1 G − 1 NOS (4311010) – M. Gilli Spring 2008 – 63 Example for solution with normal equations Least squares solution of Ax ∼ = b (normal equations solved with Cholesky factorization)     1 1 2     4 10 5 2 . 0 0 0 1 2 1     A = b = C = 10 30 11 G = 5 . 0 2 . 236 0         1 3 1     5 11 7 2 . 5 − . 670 0 . 547 1 4 1 � � 2 . 0 σ 2 = . 547 2 / 2 = 0 . 15 Solve : G ′ x = z x = − 0 . 3 � 0 . 225 � � � 0 . 500 0 − . 075 T = G − 1 = S = σ 2 T ′ T = − 1 . 118 0 . 447 − . 075 0 . 030 NOS (4311010) – M. Gilli Spring 2008 – 64 34

  34. Least squares solution with QR decomposition Length of a vector is invariant to orthogonal transformations Want to minimize � Ax − b � 2 Replace A by QR 2 QRx − b Multiply with Q ′ (orthogonal transformation) Q ′ Q Rx − Q ′ b ���� I We get � R 1 � Q ′ � � � � 2 � � � R 1 x − Q ′ 1 b � 2 2 + � Q ′ 2 b � 2 1 x − b → � � Q ′ 2 0 � � 2 2 Solve the triangular system R 1 x = Q ′ 1 b to get x The sum of squared residuals is � Q ′ 2 b � 2 2 NOS (4311010) – M. Gilli Spring 2008 – 65 Example for solution with QR decomposition A and b as before:     − . 500 . 670 . 023 . 547 � � − . 500 . 223 − . 439 − . 712 − 2 . 000 − 5 . 000     Q 1 = Q 2 = R 1 =     − . 500 − . 223 . 807 − . 217 0 − 2 . 236     − . 500 − . 670 − . 392 . 382 � − 2 . 500 � 2 . 0 � � Q ′ 1 b = c = Solve : R 1 x = c x = . 670 − . 3 � � 0 . 023 σ 2 = d ′ d/ 2 = 0 . 15 Q ′ 2 b = d = and . 547 With Matlab we simply code x = A \ b NOS (4311010) – M. Gilli Spring 2008 – 66 35

  35. Comparison between normal equations and QR In situations where � Ax − b � 2 ◦ 2 is small and κ 2 ( A ) is large the approach using QR is superior to the approach using the normal equations In situations where � Ax − b � 2 ◦ 2 is large and κ 2 ( A ) is large both approaches experience similar problems � � 2 κ 2 ( A ′ A ) = ◦ κ 2 ( A ) !!! ◦ In situations where m ≫ n the approach with normal equations uses approximatively half of elementary operations and storage compared to QR ◦ Even tough it is not possible to designate an algorithm being superior in all situations, the QR factorization is considered as the standard method (with Matlab x = A \ b ) ◦ SVD is the numerically most stable method (at the cost of a much higher operation count) NOS (4311010) – M. Gilli Spring 2008 – 67 Example for solution with SVD ? A = [1 1; 1 2; 1 3; 1 4]; b = [2 1 1 1]’; [m,n] = size(A); [U,S,V] = svd(A); c = U’*b; c1 = c(1:n); sv = diag(S); z1 = c1./sv; x = V*z1 % <-- parameters c2 = c(n+1:m); sigmac = c2 ’*c2/(m-n) % <-- sigma ^2 S = V*diag(sv .^( -2))*V’; Mcov = sigmac*S % <-- Variances and covariances for i = 1:n s(i) = V(i ,:)*V(i,:)’ * sv(i)^( -2) * sigmac; end t = x ./ sqrt(s)’ % <-- Student statistic NOS (4311010) – M. Gilli Spring 2008 – 68 36

  36. Nonlinear equations Contrarily as is the case with linear systems, the solution of nonlinear systems is a delicate matter. Generally no guarantee to converge to the true solution and the computational complexity grows very fast with the size of the system. y 1 − y 3 2 − 1 = 0 Depending on value of c there are 2 between 0 and 4 solutions!! 4 y 2 1 − y 2 + c = 0 1 1 c~0.7 c=0 0.5 0.5 y 1 y 1 0 0 −0.5 −0.5 −1 −1 −1 0 1 −1 0 1 y 2 y 2 1 1 c=−1.5 c=−1 0.5 0.5 y 1 y 1 0 0 −0.5 −0.5 −1 −1 −1 0 1 −1 0 1 y 2 y 2 NOS (4311010) – M. Gilli Spring 2008 – 69 37

  37. Nonlinear equations notation ≡ h ( y, z ) = 0 F ( y ) = 0 g 1 ( y 1 , y 2 , y 3 , y 4 , y 5 , z 2 ) = 0 g 1 ( y 1 , y 2 , y 3 , y 4 , y 5 , z 2 ) = 0 f 1 ( y 1 , y 2 , y 3 , y 4 , y 5 ) = 0 y 6 = g 2 ( y 3 , z 1 ) g 2 ( y 3 , z 1 ) − y 6 = 0 f 2 ( y 3 ) − y 6 = 0 y 3 = g 3 ( y 7 , z 4 ) g 3 ( y 7 , z 4 ) − y 3 = 0 f 3 ( y 7 ) − y 3 = 0 g 4 ( y 1 , y 2 , y 4 , y 6 , y 7 , y 8 ) = 0 g 4 ( y 1 , y 2 , y 4 , y 6 , y 7 , y 8 ) = 0 g 4 ( y 1 , y 2 , y 4 , y 6 , y 7 , y 8 ) = 0 h : g 5 ( y 5 , y 6 ) = 0 g 5 ( y 5 , y 6 ) = 0 g 5 ( y 5 , y 6 ) = 0 y 6 = g 6 ( y 7 , z 3 , z 1 ) g 6 ( y 7 , z 3 , z 1 ) − y 6 = 0 f 6 ( y 7 , z 3 ) − y 6 = 0 y 7 = g 7 ( y 3 , z 5 ) g 7 ( y 3 , z 5 ) − y 7 = 0 f 7 ( y 3 ) − y 7 = 0 g 8 ( y 3 , y 5 ) = 0 g 8 ( y 3 , y 5 ) = 0 g 8 ( y 3 , y 5 ) = 0 R n × R n → R n has derivatives in the neighborhood of the solution y ∗ ∈ R n and z ∈ R m are exogenous variables. If at least one equation is nonlinear, the system is nonlinear. In practice, often there exist a permutation of the Jacobian matrix ∂h/∂y ′ such that its structure is: In such a situation the solution of the system of equations consists in solving a sequence of recursive and interdependent systems. In the following we suppose the system to be interdependent. NOS (4311010) – M. Gilli Spring 2008 – 70 Nonlinear equations The following methods are used to solve nonlinear equations: ◦ Fixed point ◦ Newton ◦ Minimization NOS (4311010) – M. Gilli Spring 2008 – 71 Nonlinear equations (one dimension) Find value of x for which f ( x ) = 0 Also called root finding or zero finding . Algorithms: ◦ Bisection method ◦ Fixed-point iteration ◦ Newton method Compromise between slow-but-sure versus fast-but-risky . Matlab functions: fzero , roots NOS (4311010) – M. Gilli Spring 2008 – 72 Voir le cours de MN. Exercice: Find zeros of two functions one of which is a polynomial. (voir Heath) There is a compromise between slow-but-sure versus fast-but-risky . NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 72 38

  38. Systems of nonlinear equations Find vector x for which F ( x ) = 0 Algorithms: ◦ Fixed-point iteration ◦ Newton method ◦ Broyden’s method Matlab functions: fsolve NOS (4311010) – M. Gilli Spring 2008 – 73 Voir le cours de MN. Exercice: (voir Heath) NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 73 39

  39. Finding the roots of f ( x ) = 0 74 Z´ eros d’une fonction Il s’agit de trouver la valeur de x qui satisfait l’´ equation f ( x ) = 0 La solution est appel´ ee z´ ero de la fonction (il peut y avoir plusieurs solutions). Ex.: 1 − e − x = 0 . 5 Pour retrouver la formulation f ( x ) = 0 on ram` ene tous les termes ` a gauche du signe de l’´ egalit´ e 1 − e − x − 0 . 5 = 0 Dans ce cas il est facile de trouver la solution analytique x = − log(0 . 5) Souvent dans la pratique, soit il n’existe pas de solution analytique ` a ce probl` eme, soit son obtention est tr` es laborieuse → On recourt ` a la solution num´ erique NOS (4311010) – M. Gilli Spring 2008 – 74 R´ esolution graphique Soit la fonction � � x 1 = x n 1 + n > 1 et x ∈ [ x L , x U ] n − 1 6 Graphique de la fonction: f = inline(’(1 + 1/(P1-1)).^x - x.^P1’,1); xL = -2; xU = 5; x = linspace(xL,xU); plot(x,f(x,2)); 0 P1 correspond au param` etre n = 2 . On a utilis´ e les op´ erations ´ el´ ement par el´ ´ ement pour pouvoir ´ evaluer f ( x ) avec −4 x un vecteur. −0.7667 2 4 NOS (4311010) – M. Gilli Spring 2008 – 75 R´ esolution al´ eatoire Il s’agit de g´ en´ erer al´ eatoirement des valeurs pour x et de calculer la valeur de f ( x ) correspondante. Si X est l’ensemble des valeurs g´ en´ er´ ees la solution x sol est d´ efinie comme x sol = argmin | f ( x ) | x ∈ X Pour la g´ en´ eration des x on se sert de la fonction Matlab rand qui g´ en` ere des variables pseudo-al´ eatoires uniformement distribu´ ees dans l’intervalle [0 , 1[ . x L , x U et la v.a. uniforme u donn´ ees on g´ en` ere x comme x = x L + ( x U − x L ) u NOS (4311010) – M. Gilli Spring 2008 – 76 40

  40. R´ esolution al´ eatoire (cont’d) Nous avons choisi x min = − 10 et x max = 10 et avons ´ evalu´ e la fonction f ( x ) pour R = 100 000 valeurs al´ eatoires de x . Le code Matlab est le suivant: n = 2; xmin = -10; xmax = 10; R = 100000; x = xmin + (xmax - xmin)*rand(R,1); z = f(x,n); [sol,xind] = min(abs(z)); fprintf(’\n f(%8.6f) = %8.6f\n’,x(xind),sol); On a trouv´ e la solution suivante f(-0.766735) = 0.000136 (depend de la s´ equence de v.a. g´ en´ er´ ee !!) NOS (4311010) – M. Gilli Spring 2008 – 77 Recoupement par intervalles (bracketing) Il s’agit de construire des intervalles susceptibles de contenir un z´ ero. Ult´ erieurement on raffinera la recherche du z´ ero ` a l’int´ erieur de l’intervalle. On subdivise un domaine donn´ e en intervalles r´ eguliers et on examine si la fonction traverse l’axe des x en analysant le signe de la fonctions ´ evalu´ ee aux bords de l’intervalle. f ( a ) 1: ∆ = ( x max − x min ) /n 2: a = x min 3: for i = 1 : n do b = a + ∆ 4: a b if sign f ( a ) � = sign f ( b ) then 5: [ a, b ] peut contenir un z´ ero, 6: f ( b ) imprimer 7: end if a b △ a = b 8: a x min x max b 9: end for NOS (4311010) – M. Gilli Spring 2008 – 78 41

  41. Matlab code pour bracketing function ab = bracketing(f,xmin,xmax,n) k = 0; delta = (xmax - xmin)/n; a = xmin; fa = feval(f,a); for i = 1:n b = a + delta;; fb = feval(f,b); if sign(fa)~=sign(fb) k = k + 1; ab(k,:) = [a b]; end a = b; fa = fb; end NOS (4311010) – M. Gilli Spring 2008 – 79 Matlab code pour bracketing Recherche des intervalles pour la fonction g ( x ) = cos(1 /x 2 ) dans le domaine x ∈ [0 . 3 , 0 . 9] et n = 25 . f = inline(’cos(1./x.^2)’); iv = bracketing(g,0.3,0.9,25); x = linspace(0.3,0.9); subplot(211) plot(x,f(x),’k’,’LineWidth’,3), hold on for i = 1:size(iv,1) plot(iv(i,:),[0 0],’LineWidth’,20,’Color’,[1 239/255 0]) end 1 iv = 0.3000 0.3240 0.3480 0.3720 0 0.4440 0.4680 0.7800 0.8040 −1 0.2 0.4 0.6 0.8 1 NOS (4311010) – M. Gilli Spring 2008 – 80 42

  42. M´ ethode de la bisection f ( a ) z c a b f ( c ) f ( b ) Recherche le z´ ero d’une fonction dans un intervalle [ a, b ] donn´ e. On consid` ere un intervalle qui contient le z´ ero, on le divise en deux et on cherche le demi-intervalle qui contient le z´ ero. La proc´ edure est r´ eit´ er´ ee jusqu’` a ce que l’on atteint un intervalle suffisamment petit. c = a + b ≡ a + b − a (num´ eriquement plus stable) 2 2 NOS (4311010) – M. Gilli Spring 2008 – 81 Algorithme de la bisection η tol´ erance d’erreur donn´ ee erifier que f ( a ) × f ( b ) < 0 1: V´ 2: while | a − b | > η do c = a + ( b − a ) / 2 3: if f ( c ) × f ( a ) < 0 then 4: b = c ( z est ` a gauche de c ) 5: else 6: a = c ( z est ` a droite de c ) 7: end if 8: 9: end while NOS (4311010) – M. Gilli Spring 2008 – 82 Code Matlab de la bisection function c = bisection(f,a,b,tol) % Recherche zero d’une fonction avec methode de la bisection % On cherche zero dans intervalle [a, b] avec tolerance tol % La fonction doit etre de signe oppose en a et b if nargin == 3, tol = 1e-8; end fa = feval(f,a); fb = feval(f,b); if sign(fa) == sign(fb), error(’f pas de signe oppose en a et b’); end fait = 0; while abs(b-a) > 2*tol & ~fait c = a + (b - a) / 2; % Chercher centre intervalle fc = feval(f,c); % Evaluer f au centre if sign(fa) ~= sign(fc) % zero a gauche de c b = c; fb = fc; elseif sign(fc) ~= sign(fb) % zero a droite de c a = c; fa = fc; else % On tombe exactement sur zero fait = 1; end end NOS (4311010) – M. Gilli Spring 2008 – 83 43

  43. Application de la bisection g = inline(’cos(1./x.^2)’); for i = 1:size(iv,1) z(i) = bisection(g,iv(i,1),iv(i,2)); end z = 0.3016 0.3568 0.4607 0.7979 iv = 0.3000 0.3240 0.3480 0.3720 0.4440 0.4680 0.7800 0.8040 NOS (4311010) – M. Gilli Spring 2008 – 84 It´ erations de point fixe Soit f ( x ) = 0 , on isole un terme contenant x de la sorte ` a pouvoir ´ ecrire x new = g ( x old ) g est app´ ell´ ee fonction d’it´ eration L’it´ eration du point fixe est d´ efinie par l’algorithme: 1: Initialiser x (0) 2: for k = 1 , 2 , . . . jusqu’` a convergence do x ( k ) = g ( x ( k − 1) ) 3: 4: end for ◦ Notion de convergence doit ˆ etre approfondie ◦ Les valeurs successives de x ne doivent pas ˆ etre conserv´ ees dans un vecteur while ~converged x1 = g(x0) x0 = x1 end ◦ Le fait que la solution x sol v´ erifie x sol = g ( x sol ) → on appelle x sol point fixe ◦ Choix de la fonction g ( x ) est d´ eterminant pour la convergence de la m´ ethode NOS (4311010) – M. Gilli Spring 2008 – 85 44

  44. Code Matlab pour la m´ ethode du point fixe function x1 = FPI(f,x1,tol) % FPI(f,x0) Iterations du point fixe if nargin == 2, tol = 1e-8; end it = 0; itmax = 100; x0 = realmax; while ~converged1(x0,x1,tol) x0 = x1; x1 = feval(f,x0); it = it + 1; if it > itmax, error(’Maxit in FPI’); end end function rep = converged1(x0,x1,tol) rep = abs(x1-x0)/(abs(x0)+1) < tol; NOS (4311010) – M. Gilli Spring 2008 – 86 Application de la m´ ethode du point fixe Soit f ( x ) = x − x 4 / 5 − 2 = 0 et les 2 fonctions d’it´ eration g 1 ( x ) = x 4 / 5 + 2 g 2 ( x ) = ( x − 2) 5 / 4 g1 = inline(’x.^(4/5)+2’); FPI(g1,1) ans = 6.4338 g2 = inline(’(x-2).^(5/4)’); FPI(g2,8) ??? Error using ==> fpi Maxit in FPI NOS (4311010) – M. Gilli Spring 2008 – 87 45

  45. Application de la m´ ethode du point fixe (cont’d) Soit h ( x ) = x − 2 x 4 / 5 + 2 = 0 et les 2 fonctions d’it´ eration � x + 2 � 5 / 4 g 1 ( x ) = 2 x 4 / 5 − 2 g 2 ( x ) = 2 h = inline(’x - 2*x.^(4/5) + 2’); bracketing(h,0,50,30) g2 = inline(’((x+2)/2).^(5/4)’); ans = FPI(g2,18.5) 3.3333 5.0000 ans = 18.3333 20.0000 19.7603 FPI(g2,3,5) g1 = inline(’2*x.^(4/5) - 2’); ??? Error using ==> fpi FPI(g1,18.5) Maxit in FPI ans = 19.7603 NOS (4311010) – M. Gilli Spring 2008 – 88 Convergence de la m´ ethode de point fixe (illustrations) Soit f ( x ) = x − x 4 / 5 − 2 = 0 et la fonctions d’it´ eration g 1 ( x ) = x 4 / 5 + 2 3 3 3 1 1 1 5.28 4.41 3 1 3 4.41 eration g 2 ( x ) = ( x − 2) 5 / 4 Mais avec la fonctions d’it´ 9.39 8 NOS (4311010) – M. Gilli Spring 2008 – 89 46

  46. Convergence de la m´ ethode de point fixe (illustrations) eration g 1 ( x ) = 2 x 4 / 5 − 2 Soit f ( x ) = x − 2 x 4 / 5 + 2 = 0 et la fonctions d’it´ 18.64 18.5 eration g 2 ( x ) = ( x +2 2 ) 5 / 4 Mais avec la fonctions d’it´ 18.34 18.5 NOS (4311010) – M. Gilli Spring 2008 – 90 Convergence de la m´ ethode de point fixe (illustrations) Mais dans l’intervalle x ∈ [3 , 5] avec g 2 ( x ) = ( x +2 2 ) 5 / 4 4.36 3 4.5 5 et avec g 1 ( x ) = 2 x 4 / 5 − 2 on converge 4.36 3 4.5 5 NOS (4311010) – M. Gilli Spring 2008 – 91 47

  47. Convergence de la m´ ethode de point fixe erations de point fixe convergent pour x ∈ [ a, b ] si l’intervalle contient un z´ Les it´ ero et si | g ′ ( x ) | < 1 x ∈ [ a, b ] pour Si − 1 < g ′ ( x ) < 0 les it´ erations oscillent autour de la solution et si 0 < g ′ ( x ) < 1 les it´ erations convergent de fa¸ con monotone NOS (4311010) – M. Gilli Spring 2008 – 92 Convergence de la m´ ethode de point fixe (illustrations) − 1 < g ′ ( x ) < 0 g ′ ( x ) < − 1 3.92 4.37 1 1 NOS (4311010) – M. Gilli Spring 2008 – 93 Remarques: Suivant la valeur de la d´ eriv´ ee une mˆ eme fonction d’it´ eration peut converger pour certains intervalles et diverger pour d’autres !!! NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 93 M´ ethode de Newton D´ eriv´ ee ` a partir de l’expansion de la fonction f en s´ erie de Taylor: f ( x + △ x ) = f ( x )+ △ x f ′ ( x ) + R ( x ) On cherche le pas △ x pour lequel f ( x + △ x ) = 0 . Comme on neglige R ( x ) f ( x k +1 ) ≈ f ( x k ) + ( x k +1 − x k ) f ′ ( x ) = 0 f ( x k ) d’o` u f ( x k +1 ) x k +1 = x k − f ( x k ) s f ′ ( x k ) x k +1 x k � �� � s NOS (4311010) – M. Gilli Spring 2008 – 94 48

  48. Algorithme de Newton 1: Initialiser x 0 2: for k = 0 , 1 , 2 , . . . jusqu’` a convergence do x k +1 = x k − f ( x k ) 3: f ′ ( x k ) 4: end for Cr´ eer une fonction qui ´ evalue f et sa d´ eriv´ ee (2 arguments de sortie) function x1 = Newton(f,x1,tol) % Newton(f,x0) Methode de Newton if nargin == 2, tol = 1e-8; end it = 0; itmax = 10; x0 = realmax; while ~converged1(x0,x1,tol) x0 = x1; [fx,dfx] = feval(f,x0); x1 = x0 - fx/dfx; it = it + 1; if it > itmax, error(’Maxit in Newton’); end end NOS (4311010) – M. Gilli Spring 2008 – 95 Exemple d’application de l’algorithme de Newton function [f,d] = myf(x) f = exp(-x).*log(x)-x.^2+2; d = -exp(-x).*log(x)+exp(-x)./x-2*x; x1 = Newton0(’myf’,1) 1.5 1 0.5 0 −0.5 eriquement → quasi-Newton Newton avec d´ eriv´ ee approch´ e num´ 0 0.5 1 1.5 2 NOS (4311010) – M. Gilli Spring 2008 – 96 NewtonG NewtonG donne une illustration graphique de l’algorithme de Newton. Voir aussi ZeroFunctionSlides NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 96 Remarques sur la recherche du z´ ero d’une fonction ◦ Intervalle ◦ Point de d´ epart ◦ Robustesse NOS (4311010) – M. Gilli Spring 2008 – 97 49

  49. Systems of nonlinear equations 98 Iterative methods (Gauss-Seidel and Jacobi) Les m´ ethodes de Gauss-Seidel et Jacobi pour la solution de syst` emes lin´ eaires peuvent ˆ etre ´ etendus pour la r´ esolution de syst` emes non-lin´ eaires. En g´ en´ eral on normalise les ´ equations, c’est-` a-dire on les r´ e´ ecrit sous la forme y i = g i ( y 1 , . . . , y i − 1 , y i +1 , . . . , y n , z ) i = 1 , . . . , n. Rappelons que les ´ equations g i peuvent maintenant ˆ etre non-lin´ eaires. Pour la m´ ethode de Jacobi l’it´ eration g´ en´ erique s’´ ecrit alors: y ( k +1) = g i ( y ( k ) 1 , . . . , y ( k ) i − 1 , y ( k ) i +1 , . . . , y ( k ) n , z ) i = 1 , . . . , n. i a jour de y ( k +1) aussitˆ ethode de Gauss-Seidel utilise les i − 1 composantes mises ` L’it´ eration g´ en´ erique de la m´ ot qu’elles sont disponibles: y ( k +1) = g i ( y ( k +1) , . . . , y ( k +1) , y ( k ) i +1 , . . . , y ( k ) n , z ) i = 1 , . . . , n. i 1 i − 1 NOS (4311010) – M. Gilli Spring 2008 – 98 On utilise le mˆ eme crit` ere d’arrˆ et que pour le cas lin´ eaire, c’est-` a-dire on stoppera les iterations lorsque la condition suivante | y ( k +1) − y ( k ) | i i < ε i = 1 , 2 , . . . , n | y ( k ) | + 1 i est satisfaite ou ε est une tol´ erance donn´ ee. ealable car la matrice M − 1 N du th´ La convergence ne peut ˆ etre ´ etablie au pr´ eor` eme change ` a chaque it´ eration. La structure de l’algorithme it´ eratif est identique ` a celle pr´ esent´ ee pour le cas lin´ eaire: Initialiser y (1) , y (0) , ε et le nombre maximum d’it´ erations while ˜ converged( y (0) , y (1) , ε ) do y (0) = y (1) (Conserver l’it´ edente dans y (0) ) eration pr´ ec´ Evaluer y (1) avec Jacobi, Gauss-Seidel ou SOR ´ Test sur le nombre d’it´ erations end while NOS (4311010) – M. Gilli Spring 2008 – 99 50

  50. Fix point method Il s’agit d’une formulation g´ en´ erale des m´ ethodes it´ eratives ou le syst` eme d’´ equations est formalis´ e comme y = g ( y ) et l’it´ eration du point fixe s’´ ecrit y ( k +1) = g ( y ( k ) ) k = 0 , 1 , 2 , . . . epart y (0) . La condition de convergence est ´ etant donn´ e un vecteur de d´ � �   ∂g 1 ∂g 1 � � · · · � � ∂y 1 ∂y n y 1 = y sol y n = y sol 1 n � �   . . ∇ g ( y sol ) ∇ g ( y sol ) = ρ < 1  . .  . .   � � ∂g n ∂g n � � · · · � � ∂y 1 ∂y 1 y 1 = y sol y 1 = y sol 1 1 Matrice Jacobienne ´ evalu´ ee ` a la solution y sol . NOS (4311010) – M. Gilli Spring 2008 – 100 Example Illustrons la m´ ethode du point fixe en r´ esolvant le syst` eme d’´ equations non-lin´ eaires suivant: � f 1 ( y 1 , y 2 ) : y 1 + y 2 − 3 = 0 ⇔ F ( y ) = 0 f 2 ( y 1 , y 2 ) : y 2 1 + y 2 2 − 9 = 0 � 0 � 3 3 � ′ et y = 0 � ′ qu’admet ce syst` La Figure montre les deux solutions y = eme d’´ equations. 3 y 1 0 −2 0 3 y 2 � ′ et une tol´ � epart y (0) = erance de 10 − 5 , l’algorithme du point fixe En choisissant comme point de d´ 1 5 pour ce probl` eme peut s’´ ecrire comme: y1 = [1 5]’; tol = 1e-5; y0 = 1+y1; % Pour premier test dans while while ~converged(y0,y1,tol) y0 = y1; y1(1) = -y0(2) + 3; y1(2) = sqrt(-y0(1)^2 + 9); end NOS (4311010) – M. Gilli Spring 2008 – 101 51

  51. Example 0.76 0.17 y 1 −2 2.2 2.8 5 y 2 Iter Solution Erreur y ( k ) y ( k ) y ( k ) y ( k ) k 1 − y sol 2 − y sol 1 2 1 2 0 1.0000 5.0000 1.0000 2.0000 1 − 2.0000 2.8284 − 2.0000 − 0.1716 2 0.1716 2.2361 0.1716 − 0.7639 3 0.7639 2.9951 0.7639 − 0.0049 4 0.0049 2.9011 0.0049 − 0.0989 5 0.0989 3.0000 0.0989 − 0.0000 6 0.0000 2.9984 0.0000 − 0.0016 7 0.0016 3.0000 0.0016 − 0.0000 8 0.0000 3.0000 0.0000 − 0.0000 Fonctionnement de l’algorithme du point fixe avec y (0) = [ 1 5 ] function res = converged(y0,y1,tol) res = all(abs(y1-y0)./(abs(y0)+1) < tol); NOS (4311010) – M. Gilli Spring 2008 – 102 Example � ′ l’algorithme oscille entre y = � ′ et y = � � � � ′ En choisissant comme valeur initiale y = 0 2 0 0 3 3 sans converger. 3 2 y 1 0 0 3 y 2 Fonctionnement de l’algorithme du point fixe avec y (0) = [ 0 2 ] . NOS (4311010) – M. Gilli Spring 2008 – 103 52

  52. Newton method La m´ ethode de Newton pour la r´ esolution d’un syst` eme d’´ equations est une g´ en´ eralisation de l’algorithme de Newton pour la recherche du z´ ero d’une fonction unidimensionnelle. Le fait qu’un algorithme pour un probl` eme unidimensionnel puisse ˆ etre efficacement g´ en´ eralis´ e au cas multidimensionnel est une situation exceptionnelle. La m´ ethode de Newton est souvent aussi appel´ ee m´ ethode de Newton-Raphson . Dans la formulation classique on ´ ecrit le syst` eme d’´ equations comme:  f 1 ( y ) = 0   . . F ( y ) = 0 ≡ .   f n ( y ) = 0 Cette ´ ecriture est ´ equivalente ` a la notation h ( y, z ) = 0 , ` a la diff´ erence que les variables z sont directement dans F . NOS (4311010) – M. Gilli Spring 2008 – 104 Newton method e y ( k ) ∈ R n et une ´ On approche la solution y ∗ par la s´ equence { y ( k ) } k =0 , 1 , 2 ,... . ´ Etant donn´ evaluation de la matrice Jacobienne � �   ∂f 1 � ∂f 1 � · · · � � ∂y 1 ∂y n y 1 = y ( k ) y n = y ( k )   n 1  . .  ∇ F ( y ( k ) ) = . .   . .   � �   ∂f n � ∂f n � · · · � � ∂y 1 ∂y n y 1 = y ( k ) y n = y ( k ) n 1 on construit une approximation meilleure y ( k +1) de y ∗ . Pour ce faire on approche F ( y ) dans le voisinage de y ( k ) par une fonction affine: F ( y ) ≈ F ( y ( k ) ) + ∇ F ( y ( k ) )( y − y ( k ) ) . On peut r´ esoudre ce “mod` ele local” pour obtenir une valeur de y qui satisfasse F ( y ( k ) ) + ∇ F ( y ( k ) )( y − y ( k ) ) = 0 � � − 1 y = y ( k ) − ∇ F ( y ( k ) ) F ( y ( k ) ) c’est-` a-dire NOS (4311010) – M. Gilli Spring 2008 – 105 Newton method On construira alors un nouveau mod` ele local autour de la valeur obtenue pour y . A l’it´ eration k on a � � − 1 y ( k +1) = y ( k ) − ∇ F ( y ( k ) ) F ( y ( k ) ) et y ( k +1) est solution du syst` eme lin´ eaire ( y ( k +1) − y ( k ) ) ∇ F ( y ( k ) ) = − F ( y ( k ) ) � �� � � �� � � �� � s J b NOS (4311010) – M. Gilli Spring 2008 – 106 53

  53. Algorithme de Newton 1: Donner y (0) 2: for k = 0 , 1 , 2 , . . . until convergence do ´ Evaluer b = − F ( y ( k ) ) et J = ∇ F ( y ( k ) ) 3: V´ erifier la condition de J 4: J s = b solve 5: y ( k +1) = y ( k ) + s 6: 7: end for NOS (4311010) – M. Gilli Spring 2008 – 107 Newton method Illustrons la m´ ethode de Newton en r´ esolvant le syst` eme de deux ´ equations: � f 1 ( y 1 , y 2 ) : y 1 + y 2 − 3 = 0 F ( y ) = 0 ⇔ f 2 ( y 1 , y 2 ) : y 2 1 + y 2 2 − 9 = 0 et dont la matrice Jacobienne s’´ ecrit   � � ∂f 1 � ∂f 1 �   � � ∂y 1 ∂y 2 1 1 y 1 = y ( k ) y 2 = y ( k )   ∇ F ( y ( k ) ) =  . 1 2  =    � � 2 y ( k ) 2 y ( k )  ∂f 2 ∂f 2 � � 1 2 � � ∂y 1 ∂y 2 y 1 = y ( k ) y 2 = y ( k ) 1 2 � ′ alors Si l’on choisit y (0) = � 1 5 � � � � 3 1 1 F ( y (0) ) = ∇ F ( y (0) ) = et . 17 2 10 NOS (4311010) – M. Gilli Spring 2008 – 108 Newton method En r´ esolvant le syst` eme � � � � 1 1 − 3 s (0) = 2 10 − 17 � − 13 / 8 − 11 / 8 � ′ d’o` l’on obtient s (0) = u � � � � � � − . 625 0 1 1 y (1) = y (0) + s (0) = F ( y (1) ) = ∇ F ( y (1) ) = , , . 145 − 5 29 3 . 625 32 4 4 En r´ esolvant ` a nouveau le syst` eme � � � � 1 1 0 s (1) = − 5 29 − 145 4 4 32 � 145 / 272 − 145 / 272 � ′ d’o` on obtient s (1) = u � − . 092 � y (2) = y (1) + s (1) = . 3 . 092 NOS (4311010) – M. Gilli Spring 2008 – 109 54

  54. Newton method En se servant du code Matlab qui suit (en initialisant y1 ): tol = 1e-5; y0 = 1+y1; % Pour premier test dans while while ~converged(y0,y1,tol) y0 = y1; F(1) = y0(1) + y0(2) - 3; F(2) = y0(1)^2 + y0(2)^2 - 9; b = -F’; J = [1 1; 2*y0(1) 2*y0(2)]; s = J \ b; y1 = y0 + s; end on peut calculer les it´ erations jusqu’` a la convergence. NOS (4311010) – M. Gilli Spring 2008 – 110 Newton method � � ′ . Rappelons que la fonction La Figure montre comme ces it´ erations convergent vers la solution y = 0 3 Matlab converged a ´ et´ e introduite avant. 1 y 1 −0.09 −0.62 3.09 3.62 5 y 2 Iter Solution Erreur y ( k ) y ( k ) y ( k ) y ( k ) 1 − y sol 2 − y sol k 1 2 1 2 0 1.0000 5.0000 1.0000 2.0000 1 − 0.6250 3.6250 − 0.6250 0.6250 2 − 0.0919 3.0919 − 0.0919 0.0919 3 − 0.0027 3.0027 − 0.0027 0.0027 4 − 0.0000 3.0000 0.0000 0.0000 Fonctionnement de l’algorithme de Newton avec y (0) = [ 1 5 ] . NOS (4311010) – M. Gilli Spring 2008 – 111 55

  55. Newton method ethode de Newton converge aussi lorsqu’on choisit comme valeur initiale y (0) = � � ′ Remarquons que la m´ 0 2 ce qui est illustr´ e dans la Figure. 3.25 3 y 1 2 −0.25 0 y 2 Fonctionnement de l’algorithme de Newton avec y (0) = [ 0 2 ] . NOS (4311010) – M. Gilli Spring 2008 – 112 Convergence Le taux de convergence r d’une m´ ethode it´ erative est d´ efini comme � e ( k +1) � lim � e ( k ) � r = c k →∞ ou e ( k ) est l’erreur ` a l’it´ eration k et c une constante finie. On distingue alors les cas particuliers suivants: ◦ r = 1 et c < 1 , convergence lin´ eaire ◦ r > 1 , convergence super-lin´ eaire ◦ r = 2 , convergence quadratique Dans la pratique ceci correspond ` a un gain de pr´ ecision par it´ eration d’un nombre de digits constant pour une convergence lin´ eaire; l’accroissement de la pr´ ecision par it´ eration va en augmentant pour une convergence super-lin´ eaire, et dans le cas d’une convergence quadratique la pr´ ecision double ` a chaque it´ eration. NOS (4311010) – M. Gilli Spring 2008 – 113 Convergence La m´ ethode de Newton converge quadratiquement sous certaines conditions. On v´ erifie � y ( k +1) − y ∗ � βγ � y ( k ) − y ∗ � 2 ≤ k = 0 , 1 , 2 , . . . (5) e relative � ∇ F ( y ∗ ) − 1 �≤ β < 0 et γ est la constante de Lipschitz. o` u β mesure la non-lin´ earit´ La convergence est seulement garantie si y (0) se trouve dans un voisinage de y ∗ , voisinage qu’il faut d´ efinir. etriques, y (0) est d´ Pour la solution de mod` eles macro´ econom´ efini par y t − 1 ce qui constitue en g´ en´ eral un bon voisinage. La convergence quadratique des it´ erations de l’exemple peut ˆ etre v´ erifi´ ee dans le tableau de la Figure. En effet pour chaque composante y ( k ) , on v´ erifie i i | 2 . | y ( k +1) i | < | y ( k ) − y sol − y sol i i NOS (4311010) – M. Gilli Spring 2008 – 114 56

  56. Commentaire La complexit´ e de la m´ ethode de Newton est d´ etermin´ ee par la r´ esolution du syst` eme lin´ eaire. Si l’on compare les m´ ethodes it´ eratives avec la m´ ethode de Newton pour une it´ eration on remarque que la diff´ erence de travail r´ eside dans l’´ evaluation de la matrice Jacobienne et la solution d’un syst` eme lin´ eaire. Ceci peut paraˆ ıtre un d´ esavantage de la m´ ethode de Newton. Dans la pratique cependant on peut montrer que lorsque on r´ esoud le syst` eme lin´ eaire avec une m´ ethode directe creuse et lorsqu’on utilise des d´ eriv´ ees analytiques les deux m´ ethodes sont de complexit´ e comparable ´ etant donn´ e que le nombre d’it´ erations pour Newton est nettement inf´ erieur aux nombre d’it´ erations pour les m´ ethodes it´ eratives. NOS (4311010) – M. Gilli Spring 2008 – 115 Quasi-Newton La m´ ethode de Newton n´ ecessite ` a chaque it´ eration le calcul de la matrice Jacobienne, ce qui requiert evaluation de n 2 d´ eaire ( O ( n 3 ) ). Dans la situation ou l’´ l’´ eriv´ ees, et la r´ esolution d’un syst` eme lin´ evaluation des d´ eriv´ ees est coˆ uteuse on peut remplacer le calcul exact de la Jacobienne par une simple mise ` a jour. Ceci donne lieu ` a des nombreuses variantes de la m´ ethode de Newton, appellees m´ ethodes quasi-Newton . NOS (4311010) – M. Gilli Spring 2008 – 116 Broyden’s method a jour de la matrice Jacobienne se fait avec des matrices de rang un. ´ Dans ce cas la mise ` Etant donn´ e une approximation B ( k ) de la matrice Jacobienne ` a l’it´ eration k on calcule l’approximation de l’it´ eration k + 1 comme � dF ( k ) − B ( k ) s ( k ) � s ( k ) ′ B ( k +1) = B ( k ) + s ( k ) ′ s ( k ) ou dF ( k ) = F ( y ( k +1) ) − F ( y ( k ) ) et s ( k ) est la solution de B ( k ) s ( k ) = − F ( y ( k ) ) . L’algorithme de Broyden est formalis´ e ci-apr` es. NOS (4311010) – M. Gilli Spring 2008 – 117 Algorithme de Broyden 1: Donner y (0) et B (0) (approximation de ∇ F ( y (0) ) ) 2: for k = 0 , 1 , 2 , . . . until convergence do B ( k ) s ( k ) = − F ( y ( k ) ) solve 3: y ( k +1) = y ( k ) + s ( k ) 4: △ = F ( y ( k +1) ) − F ( y ( k ) ) 5: � △ − B ( k ) s ( k ) � B ( k +1) = B ( k ) + s ( k ) ′ / ( s ( k ) ′ s ( k ) ) 6: 7: end for NOS (4311010) – M. Gilli Spring 2008 – 118 57

  57. Broyden’s method Le code Matlab qui suit r´ ealise la m´ ethode de Broyden. y0 = - y1; tol = 1e-5; B = eye(2); F1 = feval(’ExSNL’,y1); while ~converged(y0,y1,tol) y0 = y1; F0 = F1; s = B \ -F0; y1 = y0 + s; F1 = feval(’ExSNL’,y1); dF = F1 - F0; B = B + ((dF - B*s)*s’)/(s’*s); end NOS (4311010) – M. Gilli Spring 2008 – 119 Broyden’s method L’´ evaluation de la fonction F ( y ) a ´ et´ e dans ce cas transf´ er´ ee dans une fonction Matlab ExSNL ce qui rend le code ind´ ependant de la r´ esolution d’un syst` eme d’´ equations particulier. function F = ExSNL(y) F = repmat(NaN,2,1); F(1) = y(1) + y(2) - 3; F(2) = y(1)^2 + y(2)^2 - 9; NOS (4311010) – M. Gilli Spring 2008 – 120 58

  58. Broyden’s method ethode de Broyden en partant avec B (0) = I et pour une valeur initiale La Figure donne les r´ esultats de m´ y (0) = [ 2 0 ] . 4.45 3 y 1 0 −1.49 0 3 5 y 2 Iter Solution Erreur y ( k ) y ( k ) y ( k ) y ( k ) 1 − y sol 2 − y sol k 1 2 1 2 0 2.0000 0 − 1.0000 0 1 3.0000 5.0000 0 5.0000 2 2.1667 0.8333 − 0.8333 0.8333 3 1.6585 1.4634 − 1.3415 1.4634 4 4.4582 − 1.4984 1.4582 − 1.4984 5 2.3304 0.6762 − 0.6696 0.6762 6 2.7384 0.2614 − 0.2616 0.2614 7 3.0847 − 0.0852 0.0847 − 0.0852 8 2.9921 0.0079 − 0.0079 0.0079 9 2.9998 0.0002 − 0.0002 0.0002 10 3.0000 0.0000 0.0000 0.0000 esultats de l’algorithme de Broyden avec B (0) = I et y (0) = [ 2 0 ] . R´ NOS (4311010) – M. Gilli Spring 2008 – 121 59

  59. Broyden’s method emarre avec B (0) = I aucune information sur la matrice Jacobienne n’est donn´ Lorsque l’algorithme d´ e et il n’est pas surprenant que la performance soit m´ ediocre (comparable ` a la m´ ethode du point fixe). Ci-apr` es on fait emarrer l’algorithme avec B (0) = ∇ F ( y (0) ) et on peut observer que dans ce cas sa performance est bien d´ sup´ erieure. Ces r´ esultats sont reproduits dans la Figure. 3 2 y 1 0 3 y 2 Iter Solution Erreur y ( k ) y ( k ) y ( k ) y ( k ) k 1 − y sol 2 − y sol 1 2 1 2 0 2.0000 0 − 1.0000 0 1 2.6250 0.3750 − 0.3750 0.3750 2 2.4044 0.5956 − 0.5956 0.5956 3 3.1100 − 0.1100 0.1100 − 0.1100 4 2.9739 0.0261 − 0.0261 0.0261 5 2.9991 0.0009 − 0.0009 0.0009 6 3.0000 0.0000 0.0000 0.0000 Algorithme de Broyden avec B (0) = ∇ F ( y (0) ) et y (0) = [ 2 0 ] . NOS (4311010) – M. Gilli Spring 2008 – 122 Damped Newton Si la valeur initiale se trouve loin de la solution, la m´ ethode de Newton et ses variantes convergent en g´ en´ eral tr` es mal. En fait la direction et surtout la longueur du pas sont peu fiables. Afin de d´ efinir un pas plus prudent on peut calculer y ( k +1) ` a l’it´ eration k comme y ( k +1) = y ( k ) + α ( k ) s ( k ) ou α ( k ) est un scalaire ` eterminer. Ainsi on choisira 0 < α ( k ) < 1 lorsqu’on est loin de la solution et α ( k ) = 1 a d´ lorsque y ( k ) est proche de la solution. Une fa¸ etre α ( k ) est de le lier ` con de contrˆ oler le param` a la valeur de � F ( y ( k ) ) � . NOS (4311010) – M. Gilli Spring 2008 – 123 60

  60. Solution by minimization Pour r´ esoudre F ( y ) = 0 on peut minimiser la fonction objectif suivante g ( y ) = � F ( y ) � p u p est une norme quelconque dans R n . Une raison qui motive cette alternative est qu’elle introduit un crit` o` ere ecider si y ( k +1) est une approximation meilleure pour y ∗ que y ( k ) . Comme ` a la solution F ( y ∗ ) = 0 on pour d´ e de comparer les vecteurs F ( y ( k +1) ) et F ( y ( k ) ) en calculant leur norme respective. Ce que l’on peut ˆ etre tent´ d´ esire est que � F ( y ( k +1) ) � p < � F ( y ( k ) ) � p ce qui conduit ` a minimiser la fonction objectif g ( y ) = 1 2 F ( y ) ′ F ( y ) min y si l’on choisit p = 2 pour la norme. On utilisera l’algorithme de Gauss-Newton ou Levenberg-Marquardt. NOS (4311010) – M. Gilli Spring 2008 – 124 Solution by minimization Graphe de � F ( y ) � p du syst` eme d’´ equations de l’exemple. NOS (4311010) – M. Gilli Spring 2008 – 125 61

  61. Classical optimization paradigm 126 Definitions X ⊂ R n ◦ constraint set (feasibility region) ◦ f : X → R objective function x ∗ ∈ X is a global minimizer of f if f ( x ∗ ) ≤ f ( x ) ∀ x ∈ X ◦ x ∗ ∈ X is a local minimizer of f if ∃ δ > 0 | f ( x ∗ ) ≤ f ( x ) ∀ x ∈ X ∩ B ( x ∗ , δ ) ◦ Minimization problem usually expressed as min x ∈ X f ( x ) Optimization is unconstraint if X = R n and constraint if X is described by equality and inequality constraints � x ∈ R n � � � c i ( x ) = 0 i ∈ E � X = c i ( x ) ≥ 0 i ∈ I NOS (4311010) – M. Gilli Spring 2008 – 126 Definitions ◦ A maximisation problem max x ∈ X f ( x ) can be reformulated into the following minimization problem min x ∈ X − f ( x ) ◦ A point x ∈ X is a feasible point ◦ A point x �∈ X is an infeasible point NOS (4311010) – M. Gilli Spring 2008 – 127 Finding an initial feasible point can be a major difficulty in optimization. NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 127 Classification Several criteria can be used for classification, e.g.: ◦ Number of variables ◦ Number of constraints ◦ Properties of the objective function: ◦ linear, quadratic, nonlinear ◦ separability, non separability ◦ sum of squares of linear or nonlinear functions ◦ convexity ◦ sparsity of the Hessian ◦ etc. NOS (4311010) – M. Gilli Spring 2008 – 128 62

  62. Classification (contd.) ◦ Properties of the set X and the constraints: ◦ convexity of X ◦ only equality or inequality constraints ◦ linear or nonlinear constraints ◦ bound constraints ℓ i ≤ x i ≤ u i ◦ sparsity of the Jacobian ◦ etc. Many other possible classifications !!! NOS (4311010) – M. Gilli Spring 2008 – 129 Many other possible classifications. NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 129 Classification (contd.) For an optimization problem with a particular structure specific algorithms have been developed taking advantage of this structure. Short list of optimization problems with special structure: ◦ Linear programming f ( x ) and c i ( x ) are linear ◦ Quadratic programming f ( x ) is quadratic and c i ( x ) are linear ◦ Convex programming f ( x ) is convex and the set X is convex � n f ( x ) = 1 j =1 f 2 j ( x ) and X = R n ◦ Nonlinear least-squares 2 ◦ Bound-constrained optimization ℓ i ≤ x i ≤ u i ◦ Network optimization: ◦ Objective function f ( x ) and constraints c i ( x ) have a special structure arising from a graph Objective function is separable f ( x ) = � n ◦ j =1 f j ( x j ) ◦ Constraints c i ( x ) are linear, involve only 1 or 2 components of x NOS (4311010) – M. Gilli Spring 2008 – 130 More definitions and notations f ( x ) , x ∈ R n ◦  ∂f  ∂x 1 .   ∇ f ( x ) = . “nabla” .   ∂f ∂x n Hessian matrix   ∂ 2 f ( x ) ∂ 2 f ( x ) ∂ 2 f ( x ) · · · ∂x 2 ∂x 1 ∂x 2 ∂x 1 ∂x n 1   ∂ 2 f ( x ) ∂ 2 f ( x ) ∂ 2 f ( x )   · · · H f ( x ) ij = ∂ 2 f ( x )  ∂x 2  ∂x 2 ∂x 1 ∂x 2 ∂x n ≡ ∇ 2 f ( x ) = 2   . . . ... ∂x i ∂x j  . . .  . . .     ∂ 2 f ( x ) ∂ 2 f ( x ) ∂ 2 f ( x ) · · · ∂x n ∂x 1 ∂x n ∂x 2 ∂x 2 n NOS (4311010) – M. Gilli Spring 2008 – 131 63

  63. Convergence rate The convergence rate of an iterative method is defined as � e ( k +1) � lim � e ( k ) � r = c k →∞ where e ( k ) the error at iteration k and c is a finite constant. We then distinguish the following situations: ◦ r = 1 et c < 1 , linear convergence ◦ r > 1 , super linear convergence ◦ r = 2 , quadratic convergence NOS (4311010) – M. Gilli Spring 2008 – 132 64

  64. Unconstraint optimization We consider the (mathematical) problem x ∈ R n f ( x ) min f : R n → R and f is assumed to be twice-continuously differentiable. First-order necessary conditions for a local minimizer x ∗ : ◦ ∇ f ( x ∗ ) = 0 Gradient (partial derivative with respect to each variable) is zero. Second-order necessary conditions for a local minimizer x ∗ : ◦ z ′ ∇ 2 f ( x ∗ ) z ≥ 0 ∀ z ∈ R n ( Hessian definite positive) f ( x ∗ + z 2 z ′ ∇ 2 f ( x ∗ + ξ ) z ) = f ( x ∗ ) + ∇ f ( x ∗ ) ′ z + 1 (Taylor expansion) � �� � � �� � x =0 NOS (4311010) – M. Gilli Spring 2008 – 133 Proof for second-order conditions with Taylor expansion. ξ ∈ [0 , 1] because we have no “Reste”. NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 133 Unconstraint optimization (contd.) Problem is related to solving nonlinear equations F ( x ) = 0 where F ( x ) corresponds to ∇ f ( x ) . Two cat´ egories of methods for solving the unconstraint optimization: ◦ Gradient-based methods ◦ Steepest descent ◦ Newton’s method ◦ Direct search methods (use only function values) ◦ Golden section method (one dimension) ◦ Simplex method (Nelder and Mead) NOS (4311010) – M. Gilli Spring 2008 – 134 Direct search methods do not mimic gradient-based methods via finite differences. NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 134 65

  65. Unconstraint optimization (One dimension) Newton’s method Consider a local quadratic approximation to the objective function (truncated Taylor series expansion) f ( x + h ) ≈ f ( x ) + f ′ ( x ) h + f ′′ ( x ) h 2 2 This quadratic function of h has its minimum for h = − f ′ ( x ) f ′′ ( x ) . This suggests the iteration scheme x ( k +1) = x ( k ) − f ′ ( x ) f ′′ ( x ) ∂f ( x + h ) /∂h = f ′ ( x ) f ′′ ( x ) + 2 f ′′ ( x ) h and the first-order condition for a minimum is We differentiate f ( x + h ) with respect to h . ∂f ( x + h ) /∂h = 0 . We recognize that this corresponds to Newton’s method for solving the nonlinear equation f ′ ( x ) = 0 (zero finding) NOS (4311010) – M. Gilli Spring 2008 – 135 We differentiate f ( x + h ) with respect to h . ∂f ( x + h ) /∂h = f ′ ( x ) f ′′ ( x ) + 2 f ′′ ( x ) h and the first-order condition for a minimum is ∂f ( x + h ) /∂h = 0 . Steepest descent has no meaning in 1 dimension. NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 135 Unconstraint optimization (One dimension) 1: Initialize x (0) ( close to the minimum !! ) 2: for k = 0 , 1 , 2 , . . . until convergence do 1.05 Compute f ′ ( x ( k ) ) and f ′′ ( x ( k ) ) 3: x ( k +1) = x ( k ) − f ′ ( x ) 4: f ′′ ( x ) 5: end for 1 0.95 0.8 1 1.2 1.4 1.6 1.8 2 2.2 x0 x0 x0 x1 x1 x0 f ( x ) = 1 − log ( x ) e − x 2 Example: f ′ ( x ) = − e − x 2 + 2 log ( x ) xe − x 2 x f ′′ ( x ) = e − x 2 + 4 e − x 2 + 2 log ( x ) e − x 2 − 4 log ( x ) x 2 e − x 2 x 2 f ( x (0) + h ) = f ( x (0) ) + f ′ ( x (0) ) h + f ′′ ( x (0) ) h 2 x1 x2 2 NOS (4311010) – M. Gilli Spring 2008 – 136 ◦ Importance of starting point !! ◦ Write Matlab code NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 136 66

  66. Unconstraint optimization (One dimension) Minimum has to be bracketed in an interval [ a, b ] → f is unimodal f ( x ) > f ( x ∗ ) if x < x ∗ f(x) f ( x ) > f ( x ∗ ) if x > x ∗ a x* b This property enables the refinement of the interval containing the solution NOS (4311010) – M. Gilli Spring 2008 – 137 NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 137 Unconstraint optimization (One dimension) √ Golden section search τ = ( 5 − 1) / 2 ≈ 0 . 618 1: Compute x 1 = a + (1 − τ )( b − a ) and f 1 = f ( x 1 ) 2: Compute x 2 = a + τ ( b − a ) and f 2 = f ( x 2 ) 3: while ( b − a ) > η do 4: if f 1 < f 2 then 5: b = x 2 6: x 2 = x 1 7: f 2 = f 1 a a a a x1 x1 x1 x2 x2 x2 b b b b 8: x 1 = a + (1 − τ )( b − a ) 9: f 1 = f ( x 1 ) 10: else 11: a = x 1 12: x 1 = x 2 a a x1 x1 x2 x2 b b 13: f 1 = f 2 14: x 2 = a + τ ( b − a ) 15: f 2 = f ( x 2 ) 16: end if 17: end while a x1 x2 b NOS (4311010) – M. Gilli Spring 2008 – 138 Golden section: A − − − − − C − − − B where CB/AB = AC/AB , (or AB 2 = BC × AC ) The particular choice of τ allows us to reduce the interval containing the minimum by a constant fraction, similar as with the bisection algorithm where this reduction is 0.5. Also we only need to compute one additional point. Exercice: Code the golden search. Application: a = 1; b = 2; x = linspace(a,b); f = inline(’1-log(x).*exp(-x.^2)’); plot(x,f(x)); hold on c = GSS(f,a,b); plot(c,f(c),’ro’) Matlab function is fminunc (analytical or numerical derivatives) NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 138 67

  67. Unconstraint optimization (multiple dimensions) Steepest decent method: −∇ f ( x ) is locally the direction of steepest descent for f The function f decreases more rapidly along the direction of the negative gradient than along any other direction Starting from an initial point x (0) the successive approximations of the solution are given by x ( k +1) = x ( k ) − α ∇ f ( x ( k ) ) where α is the solution of the one-dimensional minimization problem � � x ( k ) − α ∇ f ( x ( k ) ) min α f NOS (4311010) – M. Gilli Spring 2008 – 139 The function f decreases more rapidly along the direction of the negative gradient than along any other direction. NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 139 Unconstraint optimization (multiple dimensions) 1: Initialize x (0) 2: for k = 0 , 1 , 2 , . . . until convergence do Compute ∇ f ( x ( k ) ) 3: � � Compute α ∗ = argmin α f x ( k ) − α ∇ f ( x ( k ) ) 4: x ( k +1) = x ( k ) − α ∗ ∇ f ( x ( k ) ) 5: 6: end for Linear convergence !! NOS (4311010) – M. Gilli Spring 2008 – 140 68

  68. Unconstraint optimization (multiple dimensions) � � 2 + 0 . 05 (1 − x 1 ) 2 � � x 2 − x 2 Example: f ( x 1 , x 2 ) = exp 0 . 1 1 2 2 1 1 x 2 x 2 −1 −1 −1 1 2 −1 1 2 x 1 x 1 2 1 1.2 x 2 x ( k ) − α ∇ f ( x ( k ) ) � � f 1.04 −1 −1 1 2 1 −8 −5.09 0 x 1 α 2 1 x 2 −1 −1 1 2 x 1 NOS (4311010) – M. Gilli Spring 2008 – 141 Matlab function is fminunc (one can provide the gradient and the Hessian or approximate it numerically) NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 141 69

  69. Unconstraint optimization (multiple dimensions) Newton’s method: We consider a local quadratic approximation using a truncated Taylor series expansion: f ( x + h ) ≈ f ( x ) + ∇ f ( x ) h + 1 2 h ′ ∇ 2 f ( x ) h � �� � quadratic model First order conditions for quadratic model: ∇ f ( x ) + ∇ 2 f ( x ) h = 0 and the step h which minimizes the local model is the solution of the linear system ∇ 2 f ( x ) h = −∇ f ( x ) NOS (4311010) – M. Gilli Spring 2008 – 142 Unconstraint optimization (multiple dimensions) Newton algorithm : 1: Initialize x (0) , f ∈ C 2 2: for k = 0 , 1 , 2 , . . . until convergence do Compute ∇ f ( x ( k ) ) 3: Compute ∇ 2 f ( x ( k ) ) 4: Solve ∇ 2 f ( x ( k ) ) s ( k ) = −∇ f ( x ( k ) ) 5: x ( k +1) = x ( k ) + s ( k ) 6: 7: end for Quadratic convergence !! NOS (4311010) – M. Gilli Spring 2008 – 143 Unconstraint optimization (multiple dimensions) Quasi-Newton method : We sacrify speed of convergence in exchange of economies in computations: x ( k +1) = x ( k ) − α k H − 1 ∇ f ( x ( k ) ) k If H k = I we are in the situation of steepest descent ◦ Davidon - Fletcher - Powell ◦ Broyden - Fletcher - Goldfarb - Shannon ◦ Conjugate-gradient NOS (4311010) – M. Gilli Spring 2008 – 144 70

  70. Nonlinear least squares Minimize m g ( x ) = 1 2 r ( x ) ′ r ( x ) = 1 � r i ( x ) 2 2 i =1 ◦ r i ( x ) = y i − f ( t i , x ) i = 1 , . . . , m , (residuals) ◦ f ( t i , x ) (nonlinear model) ◦ t i (independent variable) x ∈ R n ◦ (parameters to be estimated) NOS (4311010) – M. Gilli Spring 2008 – 145 Nonlinear least squares To develop the local model we need first and second derivatives of g ( x ) m � r i ( x ) ∇ r i ( x ) = ∇ r ( x ) ′ r ( x ) ∇ g ( x ) = i =1   ∂r 1 ( x ) ∂r 1 ( x ) · · · ∂x 1 ∂x n  . .  ∇ r ( x ) = . . Jacobian matrix   . .   ∂r m ( x ) ∂r m ( x ) · · · ∂x 1 ∂x n The vector   ∂r i ( x ) ∂x 1 .   . ∇ r i ( x ) =   .   ∂r i ( x ) ∂x n corresponds to the i th row of the Jacobian matrix NOS (4311010) – M. Gilli Spring 2008 – 146 Nonlinear least squares m � � � ∇ r i ( x ) ∇ r i ( x ) ′ + r i ( x ) ∇ 2 r i ( x ) ∇ 2 g ( x ) = i =1 ∇ r ( x ) ′ ∇ r ( x ) + S ( x ) = S ( x ) = � m i =1 r i ( x ) ∇ 2 r i ( x ) NOS (4311010) – M. Gilli Spring 2008 – 147 71

  71. Example: f ( t, x ) = x 1 e x 2 t   y 1 − x 1 e x 2 t 1 Data y 2 − x 1 e x 2 t 2   t 0.0 1.0 2.0 3.0 r ( x ) =   y 3 − x 1 e x 2 t 3   y 2.0 0.7 0.3 0.1 y 4 − x 1 e x 2 t 4 Jacobian     ∂r 1 ( x ) ∂r 1 ( x ) − e x 2 t 1 − x 1 t 1 e x 2 t 1 ∂x 1 ∂x 2     ∂r 2 ( x ) ∂r 2 ( x ) − e x 2 t 2 − x 1 t 2 e x 2 t 2     ∂x 1 ∂x 2     ∇ r ( x ) = =     ∂r 3 ( x ) ∂r 3 ( x ) − e x 2 t 3 − x 1 t 3 e x 2 t 3     ∂x 1 ∂x 2     ∂r 4 ( x ) ∂r 4 ( x ) − e x 2 t 4 − x 1 t 4 e x 2 t 4 ∂x 1 ∂x 2 NOS (4311010) – M. Gilli Spring 2008 – 148 Example: f ( t, x ) = x 1 e x 2 t First order derivatives      − � 4 ∂g ( x ) i =1 r i ( x ) e x 2 t i ∂x 1  = ∇ r ( x ) ′ r ( x ) = ∇ g ( x ) =   − � 4 ∂g ( x ) i =1 r i ( x ) x 1 t i e x 2 t i ∂x 2 m second order derivatives     ∂ 2 r i ( x ) ∂ 2 r i ( x ) − t i e x 2 t i  0 ∂x 2 ∂x 1 ∂x 2 ∇ 2 r i ( x ) =  = 1   ∂ 2 r i ( x ) ∂ 2 r i ( x ) − t i e x 2 t i − x 1 t 2 i e x 2 t i ∂x 2 ∂x 1 ∂x 2 2   4 4 � � ( e x 2 t i ) 2 x 1 t i ( e x 2 t i ) 2     4 t i e x 2 t i  0 �    . ∇ 2 g ( x ) = i =1 i =1 − r i ( x )   4 4   t i e x 2 t i x 1 t 2 i e x 2 t i � � i =1  x 1 ( t i e x 2 t i ) 2 ( x 1 t i e x 2 t i ) 2  i =1 i =1 NOS (4311010) – M. Gilli Spring 2008 – 149 72

  72. Quadratic approximation m c ( x ) Approximate g ( x ) in the neighborhood of x c m c ( x ) = g ( x c ) + ∇ g ( x c ) ′ ( x − x c ) + 1 2 ( x − x c ) ′ ∇ 2 g ( x c )( x − x c ) Seek x + = x c + s N satisfying first order conditions ( ∇ m c ( x + ) = 0 ) ∇ m c ( x + ) = ∇ g ( x c ) + ∇ 2 g ( x c ) ( x + − x c ) = 0 � �� � s N s N is the Newton path. Solution of linear system ∇ 2 g ( x c ) s N = −∇ g ( x c ) 1: Initialize x (0) 2: for k = 0 , 1 , 2 , . . . until convergence do ∇ 2 g ( x ( k ) ) s ( k ) = −∇ g ( x ( k ) ) 3: Solve N x ( k +1) = x ( k ) + s ( k ) 4: N 5: end for NOS (4311010) – M. Gilli Spring 2008 – 150 Detail: When deriving m c ( x ) with respect to x , g ( x c ) vanishes, ∇ g ( x c ) ′ ( x − x c ) becomes ∇ g ( x c ) because ∂Ax ∂x = A ′ 2 ( x − x c ) ′ ∇ 2 g ( x c )( x − x c ) becomes ∇ 2 g ( x c ) ( x + − x c ) because ∂x ′ Ax and = 2 Ax 1 ∂x NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 150 Approximations for ∇ 2 g ( x ) S ( x ) = � m ∇ 2 g ( x ) = ∇ r ( x ) ′ ∇ r ( x ) + S ( x ) i =1 r i ( x ) ∇ 2 r i ( x ) Gauss-Newton method (ignore S ( x ) ) 1: Initialize x (0) 2: for k = 0 , 1 , 2 , . . . until convergence do Compute r ( x ( k ) ) and ∇ r ( x ( k ) ) 3: � � ∇ r ( x ( k ) ) ′ ∇ r ( x ( k ) ) s ( k ) GN = −∇ r ( x ( k ) ) ′ r ( x ( k ) ) Solve 4: x ( k +1) = x ( k ) + s ( k ) 5: GN 6: end for Normal equations in statement 4 can be considered as the overidentified system: ∇ r ( x ( k ) ) s ( k ) GN ≈ − r ( x ( k ) ) (Solved with QR) Problem if ∇ r ( x ( k ) ) has not full rank !!! NOS (4311010) – M. Gilli Spring 2008 – 151 73

  73. Approximations for ∇ 2 g ( x ) Levenberg-Marquardt method (replace S ( x ) by µI ) 1: Initialize x (0) 2: for k = 0 , 1 , 2 , . . . until convergence do Compute r ( x ( k ) ) and ∇ r ( x ( k ) ) 3: � � ∇ r ( x ( k ) ) ′ ∇ r ( x ( k ) ) + µI s ( k ) LM = −∇ r ( x ( k ) ) ′ r ( x ( k ) ) Solve 4: x ( k +1) = x ( k ) + s ( k ) 5: LM 6: end for Normal equations in statement 4 A ′ Ax A ′ b = � � � � � � ∇ r ( x ( k ) ) � � r ( x ( k ) ) 1 s ( k ) 1 ∇ r ( x ( k ) ) ′ = ∇ r ( x ( k ) ) ′ µ 2 I − µ 2 I 1 LM 0 2 I µ Solved as the overidentified system: � � � � ∇ r ( x ( k ) ) r ( x ( k ) ) s ( k ) LM ≈ − 1 0 µ 2 I Choice for 0 ≤ µ < ∞ : In practice ≈ 10 − 2 ◦ µ = 0 → Gauss-Newton ◦ µ large → steepest descent ◦ Levenberg-Marquardt is very robust !! NOS (4311010) – M. Gilli Spring 2008 – 152 74

  74. 75

  75. Example: f ( t, x ) = x 1 e x 2 t f(t,x) = x1 * exp( x2 * t ) f(t,x) = x1 * exp( x2 * t ) 2 2 x1=2.5 x2=−2.5 y y 0.7 0.7 0.3 0.3 0.1 0.1 −1 0 1 2 3 4 −1 0 1 2 3 4 t t f(t,x) = x1 * exp( x2 * t ) f(t,x) = x1 * exp( x2 * t ) 2 x1=2.5 x2=−2.5 2 x1=2.5 x2=−2.5 x1=2.1 x2=−2.9 y y 0.7 0.7 0.3 0.3 0.1 0.1 −1 0 1 2 3 4 −1 0 1 2 3 4 t t f(t,x) = x1 * exp( x2 * t ) 2 x1=2.5 x2=−2.5 2 x1=2.1 x2=−2.9 1.5 y 1 0.7 0.3 0.5 0.1 0 0 4 −2 2 −1 0 1 2 3 4 t −4 0 0 0 −0.5 −0.5 −1 −1 −1.5 −1.5 −2 −2 −2.5 −2.5 −3 −3 −3.5 −3.5 −4 −4 0 1 2 3 4 0 1 2 3 4 0 76 −0.5 −1 −1.5 −2

  76. mcnldemo.m Use ”contour”, ”good” and [1.6 -2.2] NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 153 77

  77. Direct search methods In order to be successful (rapid and global convergence), gradient based methods, require f to be “ nice ” (practice is not always so rosy !!) We consider problems with the following feature: ◦ Derivatives of f ( x ) are not available, or do not exist ◦ Computation of f ( x ) is very expensive (e.g. obtained by a huge simulation) ◦ Values of f are inherently inexact or “ noisy ” (affected by Monte Carlo variance) ◦ We are only interested in an improvement of f rather than in a fully accurate optimum This class of problems can be approached with direct search methods NOS (4311010) – M. Gilli Spring 2008 – 154 NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 154 Direct search methods Historical note: ◦ First suggested in the 1950s ◦ Development until mid 1960s ((Hooke and Jeeves, 1961), (Spendley et al. , 1962), (Nelder and Mead, 1965)), considered as part of mainstream optimization ◦ By the 1980s direct search methods became invisible in the optimization community (However remained extremely popular among practitioners, in particular in chemistry, chemical engineering, medicine, etc.) ◦ Regained interest due to the work by (Torczon, 1989) ◦ Actual research by (Torczon, 1997), (Lagarias et al. , 1999), (Frimannslund and Steihaug, 2004), etc. NOS (4311010) – M. Gilli Spring 2008 – 155 78

  78. Direct search methods Simplex-based direct search introduced by Spendley et al. (1962): f is evaluated at the vertices of a simplex in the search space. The simplex evolves toward the minimum by constructing a new simplex obtained by reflecting the original simplex away from the vertex with the largest value of f . NOS (4311010) – M. Gilli Spring 2008 – 156 Nelder-Mead algorithm   0 d 1 d 2 · · · d 2 0 d 2 d 1 · · · d 2   x (3) X   n × ( n +1) = . . . . ... . . . .   . . . .   x (2) s 0 d 2 d 2 · · · d 1 √ d 1 = s ( √ n + 1 + n − 1) / ( n x (1) x ¯ 2) √ d 2 = s ( √ n + 1 − 1) / ( n 2) x R + NOS (4311010) – M. Gilli Spring 2008 – 157 79

  79. Nelder-Mead algorithm At iteration k simplex is defined by vertices x ( i ) , i = 1 , . . . , n + 1 ◦ Vertices are renamed such that f ( x (1) ) ≤ f ( x (2) ) ≤ . . . ≤ f ( x ( n +1) ) ◦ � n x = 1 i =1 x ( i ) ◦ Compute mean of all vertices except the worst ¯ i = 1 , . . . , n n f ( x (1) ) f ( x (3) ) f ( x (3) ) f ( x (2) ) f ( x (2) ) f ( x (1) ) x (3) x (2) ¯ x x (1) x (3) x (2) x (1) NOS (4311010) – M. Gilli Spring 2008 – 158 Nelder-Mead algorithm (cont’d) x R = (1 + ρ ) ¯ x − ρ x ( n +1) ◦ Compute reflection x E = (1 + ρ ) x R − ρ ¯ If f ( x ( R) ) < f ( x (1) ) compute expansion ◦ x f ( x (3) ) f ( x (3) ) f ( x (2) ) f ( x (2) ) f ( x R ) f ( x R ) f ( x (1) ) f ( x (1) ) x (2) x (2) x ¯ x ¯ x (3) x (3) x R x R x E x (1) x (1) NOS (4311010) – M. Gilli Spring 2008 – 159 When computing the expansion we do not check whether the function goes up again! NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 159 80

  80. Nelder-Mead algorithm (cont’d) If f ( x ( R) ) < f ( x ( n +1) ) but f ( x ( R) ) ≥ f ( x ( n ) ) ◦ x O = (1 + ψρ ) ¯ x − ψρ x ( n +1) compute out-contraction f ( x (3) ) f ( x (3) ) f ( x (2) ) f ( x R ) f ( x (2) ) f ( x R ) f ( x (1) ) f ( x (1) ) x (2) x (2) x O x ¯ ¯ x x (3) x (3) x R x R x (1) x (1) NOS (4311010) – M. Gilli Spring 2008 – 160 Nelder-Mead algorithm (cont’d) If f ( x ( R) ) > f ( x ( n +1) ) ◦ x I = (1 − ψρ ) ¯ x + ψρ x ( n +1) compute in-contraction f ( x R ) f ( x R ) f ( x (3) ) f ( x (3) ) f ( x (2) ) f ( x (2) ) f ( x (1) ) f ( x (1) ) x (2) x (2) x I x ¯ ¯ x x (3) x (3) x R x R x (1) x (1) NOS (4311010) – M. Gilli Spring 2008 – 161 81

  81. Nelder-Mead algorithm (cont’d) If outside or inside contraction results in no improvement over f ( x ( n +1) ) ◦ x ( i ) = x (1) − σ ( x ( i ) − x (1) ) simplex shrinks i = 2 , . . . , n + 1 f ( x I ) f ( x O ) f ( x (3) ) f ( x (1) ) x (2) x I x (2) x O x ¯ x (3) x R x (3) x (1) x (1) NOS (4311010) – M. Gilli Spring 2008 – 162 Nelder-Mead algorithm (cont’d) 1: Construct vertices x (1) , . . . , x ( n +1) of starting simplex 2: repeat Rename vertices such that f ( x (1) ) ≤ . . . ≤ f ( x ( n +1) ) 3: if f ( x R ) < f ( x (1) ) then 4: if f ( x E ) < f ( x R ) then x ∗ = x E else x ∗ = x R 5: 6: else if f ( x R ) < f ( x ( n ) ) then 7: x ∗ = x R 8: 9: else if f ( x R ) < f ( x ( n +1) ) then 10: if f ( x O ) < f ( x ( n +1) ) then x ∗ = x O else shrink 11: 12: else if f ( x I ) < f ( x ( n +1) ) then x ∗ = x I else shrink 13: 14: end if 15: end if 16: end if if not shrink then x ( n +1) = x ∗ (Replace worst vertex by x ∗ ) 17: 18: until stopping criteria verified NOS (4311010) – M. Gilli Spring 2008 – 163 82

  82. Nelder-Mead algorithm (cont’d) Matlab function fminsearch implements the Nelder-Mead algorithm. Example: Minimize the function 1 + x 2 − 11) 2 + ( x 1 + x 2 f ( x ) = ( x 2 2 − 7) 2 and choose x = [0 − 2] ′ as starting point. F = inline (’(x(1)^2 + x(2) -11)^2 + (x(1) + x(2)^2 - 7)^2 ’); x0 = [0 -2]; options = optimset(’fminsearch ’); [x,f,FLAG ,output] = fminsearch(F,x0 ,options ); If F is defined by a function with arguments a1,a2,... then the call is fminsearch(@F,x0,options,a1,a2,...) Similar algorithms are multidirectional search and pattern search. NOS (4311010) – M. Gilli Spring 2008 – 164 More examples Give complete example where the objective function has to be computed with a function NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 164 Examples Matlab demo with DirectSearchDemo.m : ◦ prob283 ◦ Tankership size NOS (4311010) – M. Gilli Spring 2008 – 165 Cost function for refined oil ($/kl) c = cc Crude oil price ($/kl) + ci Insurance cost ($/kl) + cx Droit de douane ($/kl) 2 . 09 × 10 4 t − 0 . 3017 � � + / 360 Freit + 1 . 064 × 10 6 a t 0 . 4925 Chargement et d´ echargement 52 . 47 q × 360 + 4 . 242 × 10 4 a t 0 . 7952 + 1 . 81 i p ( n t + 1 . 2 q ) Frais de mouillage 52 . 47 q × 360 + 4 . 25 × 10 3 a ( n t + 1 . 2 q ) Submarine pipe cost 52 . 47 q × 360 + , 5 . 042 × 10 3 q − 0 . 1899 Tank area cost 360 + 0 . 1049 q 0 . 671 Refining 360 NOS (4311010) – M. Gilli Spring 2008 – 166 83

  83. Cost function for refined oil ($/kl) (cont’d) a 0 . 2 Charges fixes annuelles (en pourcent) cc 12 . 5 Prix du brut ($/kl) ci 0 . 5 Coˆ ut d’assurance ($/kl) cx 0 . 9 Droits de douane ($/kl) i . 10 Taux d’int´ erˆ et n 2 Nombre de ports Prix du terrain ($/m 2 ) p . 7 q Capacit´ e de raffinage (bbl/jour) (1 kl = 6.29 bbl) t Tonnage du p´ etrolier (kl) c Coˆ ut du petrol raffin´ e ($/kl) Solution: q = 174 835 t = 484 687 NOS (4311010) – M. Gilli Spring 2008 – 167 84

  84. Optimization heuristics 168 Outline ◦ Standard optimization paradigm Limits ◦ Heuristic optimization paradigm ◦ Heuristics ◦ Stochastics of the solution ◦ Classification Examples NOS (4311010) – M. Gilli Spring 2008 – 168 Motivation ◦ In many fields the scope of quantitative research broadens due to increasing availability of data sets (high frequency time series, large panels, etc.) ◦ To extract relevant information from these data new statistical procedures are developed ◦ Often these procedures result in highly complex optimization problems which cannot be tackled by standard algorithms (existence of many local optima, non existent derivatives, discontinuities, etc.) ◦ Most applications circumvent these problems by simplifying models until they fit the traditional optimization paradigm → we loose relevant aspects or using ad hoc procedures → quality of result is uncertain NOS (4311010) – M. Gilli Spring 2008 – 169 Motivation (contd.) ◦ Need for further dissemination of recent advances in heuristic optimization tools (suitable to deal with highly complex optimization problems) ◦ Use of such tools is in many fields still limited ◦ Neglect due to: missing information about computational complexity of the relevant optimization problems, missing knowledge about optimization heuristics, missing formal framework for the analysis of the additional stochastic component introduced by the heuristics NOS (4311010) – M. Gilli Spring 2008 – 170 85

  85. Motivation (contd.) Promoting: ◦ use of optimization heuristics ◦ standardizing the formal framework Should contribute that a solution (to a complex non convex optimization problem), solved with a heuristic, correctly described, gets the same consideration as the solution of a convex problem obtained with a standard method. Understand the random character of the solution of a heuristic method (similar as with results from estimation) NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 171 NOS (4311010) – M. Gilli Spring 2008 – 171 The standard optimization paradigm and its limits Optimization problems in estimation and modeling typically expressed as: max x ∈X f ( x ) (6) (6) often synonymous with solution x opt assumed to exist and frequently to be unique !! McCullough and Vinod (1999, p. 635) state: ‘ Many textbooks convey the impression that all one has to do is use a computer to solve the problem, the implicit and unwarranted assumptions being that the computer’s solution is accurate and that one software package is as good as any other ’. NOS (4311010) – M. Gilli Spring 2008 – 172 Globally convex problems 18 17.95 17.9 17.85 40 17.8 17.75 20 0.2 0.705 0.71 0.15 0.695 0.7 0 0.1 0.69 −3 −2 −1 0 1 2 3 First order conditions provide the solution !! NOS (4311010) – M. Gilli Spring 2008 – 173 86

  86. Locally convex problems 1.05 1 x0 0.95 0.8 1 1.2 1.4 1.6 1.8 2 2.2 x0 x1 x0 2 1.5 1 0.5 0 0 4 −2 2 −4 0 NOS (4311010) – M. Gilli Spring 2008 – 174 mcnldemo.m Use ”contour”, ”good” and [1.6 -2.2] NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 174 Problems with many local minima !! Where do such problems arise ?? NOS (4311010) – M. Gilli Spring 2008 – 175 87

  87. Robust LSQ 8 8 6 6 4 4 2 2 0 0 0 2 4 6 8 0 2 4 6 8 0.22 0.2 0.18 β 0 0.16 8 0.14 6 0.12 4 0.1 2 0.69 0.695 0.7 0.705 0.71 β 1 0 0 2 4 6 8 Least median of squares estimator (LMS) y i = x T i θ + ǫ i i = 1 , . . . , N ˆ θ LMS = argmin Q LMS ( θ ) θ Q LMS ( θ ) = med i ( r 2 i ) median of squared residuals r 2 i = ( y i − x T i θ ) 2 300 6 200 4 100 2 0 100 2 10 0 1 5 50 2 0 0 1 0 −5 −1 −1 θ 1 0 θ 2 −10 θ 1 θ 2 NOS (4311010) – M. Gilli Spring 2008 – 176 88

  88. Example from finance VaR minimization Minimize (1 − β ) -quantile of the empirical distribution of losses of portfolio 1 1−beta 0 VaR NOS (4311010) – M. Gilli Spring 2008 – 177 Objective function for VaR minimization 5 x 10 4.2 4 3.8 3.6 3.4 4 3.5 3.2 3 4 12 x 10 11 2.5 10 9 8 2 4 7 x 10 x 1 x 2 NOS (4311010) – M. Gilli Spring 2008 – 178 Note on portfolio examples Example with 3 assets. The third asset is determined by the budget constraint. The number of assets is expressed as integer variables. Only objective functions for two dimensional problems can be illustrated. In real applications it is most likely that we have to optimize with respect to many variables, which makes the problem much more complex. NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 178 89

  89. Ω minimization 1 I 2 p Ω = I 2 I 1 I 1 0 0 NOS (4311010) – M. Gilli Spring 2008 – 179 Objective function for Ω minimization 0.36 0.34 0.32 0.3 0.28 0.26 5 4 0.6 0.8 3 1 1.2 1.4 2 4 x 1 5 x 2 x 10 x 10 NOS (4311010) – M. Gilli Spring 2008 – 180 Model selection Discrete problem: Given K possible explanatory variables select the subset which best fits the observations (minimizes a given objective function) search space Ω = { 0 , 1 } K NOS (4311010) – M. Gilli Spring 2008 – 181 90

  90. Model selection Mixed (discrete and continuous) problem: VAR(p) model p = 2   µ 1 · · · µ 3 φ 1 11 · · · φ 1   n 1     y 1 · · · y n 1 y 1 · · · y n y 1 · · · y n     . .   3 3 2 2 1 1 . .   . .   . . . . .  . .      . . . . .  . .    φ 1 1 n · · · φ 1 . . . . .   . .     =  nn         y 1 T − 1 · · · y n   1 y 1 T − 2 · · · y n T − 2 y 1 T − 3 · · · y n  φ 2 11 · · · φ 2      T − 3  T − 1 n 1     y 1 · · · y n 1 y 1 T − 1 · · · y n T − 1 y 1 T − 2 · · · y n . .   T − 2 T T . .   . .     φ 2 1 n · · · φ 2 nn ℓ lag φ ℓ Y = X B + E i lhs var. ij j rhs var. ˆ B = ( X ′ X ) − 1 X ′ Y estimate B : (OLS) NOS (4311010) – M. Gilli Spring 2008 – 182 Model selection Estimation in mixture models: Mixture of two (normal) distributions with proportion p y i = z i u 1 i + (1 − z i ) u 2 i i = 1 , . . . , n z i ∼ B (1 , p ) are binomial and u 1 ∼ N ( µ 1 , σ 2 1 ) et u 2 ∼ N ( µ 2 , σ 2 2 ) � � � � � � (1 − p ) f 2 ( y i ; µ 2 , σ 2 ) LL = log p f 1 ( y i ; µ 1 , σ 1 ) + log i ∈{ i | z i =1 } i ∈{ i | z i =0 } f j ( y i ; µ j , σ j ) , j = 1 , 2 is the normal density for y i Minimization with respect to µ 1 , σ 1 , µ 2 , σ 2 and p NOS (4311010) – M. Gilli Spring 2008 – 183 91

  91. Model selection 9000 8900 8800 8700 8600 8500 50.5 60 50 40 20 49.5 µ 1 0 z NOS (4311010) – M. Gilli Spring 2008 – 184 Minimization of functions generated by Monte Carlo simulations Estimation of parameters of an agent based model of financial market 4 3 2 1 0 0.5 0.4 0.01 0.02 0.03 0.04 0.05 0.3 0.2 0.1 δ ε NOS (4311010) – M. Gilli Spring 2008 – 185 Cluster analysis Objective function with all sort of criteria. NOS (4311010) – M. Gilli Spring 2008 – 186 92

  92. Limits of the classical optimization paradigm Classical optimization paradigm understood as: ◦ Solution identified by means of enumeration or differential calculus ◦ Existence of (unique) solution presumed ◦ Convergence of classical optimization methods for the solution of the corresponding first-order conditions Many optimization problems in statistics fall within this category (e.g. OLS estimation) However many optimization problems resist this standard approach Limits of the classical optimization: ◦ Problems which do not fulfill the requirements of these methods. ◦ Cases where the standard optimization paradigm can be applied, but problem sizes may hinder efficient calculation. NOS (4311010) – M. Gilli Spring 2008 – 187 Classification of estimation and modelling problems (relative to the classical optimization paradigm) Universe of estimation and modelling problems: Discrete Continuous NOS (4311010) – M. Gilli Spring 2008 – 188 Universe of problems Gray: Set X of possible solutions: ◦ continuous ◦ discrete Green: Easy to solve ◦ continuous: (e.g. LS estimation) allowing for an analytical solution ◦ discrete: allowing for a solution by enumeration for small scaled problems Red: Tractable by standard approximation methods: Solution can be approximated reasonably well by standard algorithms (e.g. gradient methods) Complementary set: Straightforward application of standard methods will, in general, not even provide a good approximation of the global optimum NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 188 93

  93. The heuristic optimization paradigm Methods: ◦ Based on concepts found in nature ◦ Have become feasible as a consequence of growing computational power ◦ Although aiming at high quality solution, they cannot pretend to produce the exact solution in every case with certainty Nevertheless , a stochastic high–quality approximation of a global optimum is probably more valuable than a deterministic poor–quality local minimum provided by a classical method or no solution at all. ◦ Easy to implement to different problems ◦ Side constraints on the solution can be taken into account at low additional cost NOS (4311010) – M. Gilli Spring 2008 – 189 Overview of optimization heuristics Two broad classes: ◦ Construction methods (greedy algorithms) 1: A = { a 1 , a 2 , . . . , a n } components of solution 2: E = ∅ (solution), S (set of feasible solutions) 3: while solution not complete do Choose a ∈ A (following a given ordre) 4: if E ∪ { a } ∈ S then 5: E = E ∪ { a } 6: A = A \{ a } 7: end if 8: 9: end while Ex.: Construction of maximum spanning tree ◦ Local search methods Solution space not explored systematically A particular heuristic is characterized by the way the walk through the solution domain is organized NOS (4311010) – M. Gilli Spring 2008 – 190 Submodular functions The greedy algorithm finds a global optimum for submodular functions. Let N be a finite set (the ground set 2 N denotes the set of all subsets of N A set function f : 2 N → R is submodular iff for all A, B ⊆ N f ( A ∪ B ) + f ( A ∩ B ) ≤ f ( A ) + f ( B ) stmv@adk.commerce.ubc.ca NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 190 94

  94. Classical local search for minimization 1: Generate current solution x c 2: while stopping criteria not met do Select x n ∈ N ( x c ) 3: (neighbor to current sol.) if f ( x n ) < f ( x c ) then x c = x n 4: 5: end while Selection of neighbor x n and criteria for acceptance define different walks through the solution space Stopping criteria (a given number of iterations) NOS (4311010) – M. Gilli Spring 2008 – 191 Classical meta-heuristics • Trajectory methods ◦ Threshold methods (TM) − Simulated annealing (SA) − Threshold accepting (TA) ◦ Tabu search (TS) • Population based methods ◦ Genetic algorithms (GA) ◦ Ant colonies (AC) ◦ Differential evolution (DE) NOS (4311010) – M. Gilli Spring 2008 – 192 Classical meta-heuristics Different rules for choice and/or acceptance of neighbor solution All accept uphill moves (in order to escape local minima) 0.3519 NOS (4311010) – M. Gilli Spring 2008 – 193 0.249 0.0683 0.2768 Illustration of neighborhood choice with Matlab script C:/Projects/Heuropt/CyprusLecture/SearchEx00.m NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 193 Trajectory methods Simulated annealing (SA) ◦ Kirkpatrick et al. (1983) ◦ Imitates the annealing process of solids Improvement of solution for move from x c to x n always accepted ◦ ◦ Accepts uphill move, only with given probability (decreases in a number of rounds to zero) NOS (4311010) – M. Gilli Spring 2008 – 194 95

  95. SA (contd.) 1: Generate current solution x c , initialize R and T 2: for r = 1 to R do while stopping criteria not met do 3: Compute x n ∈ N ( x c ) (neighbor to current sol.) 4: Compute △ = f ( x n ) − f ( x c ) and generate u (urv) 5: ( e − △ /T > u ) then x c = x n if ( △ < 0) or 6: end while 7: e − ∆ /T e − ∆ /T Reduce T 8: 9: end for 1 1 0.5 0 1 u 1 0.5 0.5 ∆ /T 0 0 T ∆ NOS (4311010) – M. Gilli Spring 2008 – 195 Threshold accepting (TA) Dueck and Scheuer (1990) Deterministic analog of Simulated Annealing ◦ Sequence of temperatures T replaced by sequence of thresholds τ . ◦ Statement 6. of SA algorithm becomes: then x c = x n if △ < τ ◦ Statement 8: threshold τ reduced instead of T 6: Compute △ = f ( x n ) − f ( x c ) and generate u (urv) ( e − △ /T > u ) then x c = x n if △ < τ then x c = x n 7: if ( △ < 0) or 8: end while 9: Reduce T Reduce τ NOS (4311010) – M. Gilli Spring 2008 – 196 Tabu search (TS) ◦ Glover and Laguna (1997) ◦ Designed for exploration of discrete search spaces with finite set of neighbor solutions ◦ Avoids cycling (visiting same solution more than once) by use of short term memory (tabu list, most recently visited solutions) NOS (4311010) – M. Gilli Spring 2008 – 197 96

  96. TS (contd.) 1: Generate current solution x c and initialize tabu list T = ∅ 2: while stopping criteria not met do Compute V = { x | x ∈ N ( x c ) }\ T 3: Select x n = min( V ) 4: x c = x n and T = T ∪ x n 5: Update memory 6: 7: end while 2: Stopping criterion: given number of iterations or number of consecutive iterations without improvement Choice of x n may or may not investigate all neighbors of x c 3: If more than one element is considered, x n corresponds to the best neighbor solution 5: A simple way to update memory is to remove older entries from tabu list (circular memory) NOS (4311010) – M. Gilli Spring 2008 – 198 Population based methods Genetic algorithm (GA) ◦ Holland (1975) ◦ Imitates evolutionary process of species that sexually reproduce. ◦ Do not operate on a single current solution, but on a set of current solutions ( population ). New individuals P ′′ generated with cross-over : ◦ combines part of genetic patrimony of each parent and applies a random mutation If new individual ( child ), inherits good characteristics from parents → higher probability to survive. NOS (4311010) – M. Gilli Spring 2008 – 199 GA (contd.) 1: Generate current population P of solutions 2: while stopping criteria not met do Select P ′ ⊂ P (mating pool), initialize P ′′ = ∅ (childs) 3: 4: for i = 1 to n do Select individuals x a and x b at random from P ′ 5: Apply cross-over to x a and x b to produce x child 6: Randomly mutate produced child x child 7: P ′′ = P ′′ ∪ x child x a : 1 0 1 0 1 1 1 0 1 0 1 1 0 1 1 1 0 1 0 1 1 1 1 1 0 8: end for 9: x child : 1 0 1 0 1 1 1 0 1 0 1 1 0 1 1 1 0 0 0 1 0 1 0 0 1 P = survive( P ′ , P ′′ ) 0 10: 11: end while x b : 0 0 1 1 0 1 1 1 0 0 1 1 0 0 1 0 0 0 0 1 0 1 0 0 1 3: Set of starting solutions 4–10: Construction of neighbor solution NOS (4311010) – M. Gilli Spring 2008 – 200 97

  97. GA (contd.) Select P ′ ⊂ P (mating pool), initialize P ′′ = ∅ (childs) for i = 1 to n do Compute P ′′ end for P = survive( P ′ , P ′′ ) Survivors P (new population) formed either by: last generated individuals P ′′ (childs) ◦ P ′′ ∪ { fittest from P ′ } ◦ only the fittest from P ′′ ◦ the fittest from P ′ ∪ P ′′ ◦ NOS (4311010) – M. Gilli Spring 2008 – 201 Ant systems, Ant Colony Optimization (AS, ACO) Colorni et al. (1992a) ◦ Imitates way ants search food and find their way back to their nest. ◦ First an ant explores its neighborhood randomly. As soon as a source of food is found it starts to transport food to the nest leaving traces of pheromone on the ground. The traces guide other ants to the source. ◦ Intensity of the pheromone traces depend on quantity and quality of food available at source and from distance between source and nest. As for a short distance more ants will travel on the same trail in a given time interval the process is reinforced. ◦ Ants preferably travel along important trails. → Their behavior is able to optimize their work. ◦ Pheromone trails evaporate. Once a source of food is exhausted the trails will disappear and the ants will start to search for other sources. NOS (4311010) – M. Gilli Spring 2008 – 202 98

  98. AS (contd.) F N NOS (4311010) – M. Gilli Spring 2008 – 203 AS (contd.) Reinforced process: On shorter route more ants can pass within same time ⇓ more pheromone on shorter route ⇓ ⇑ more ants attracted NOS (4311010) – M. Gilli Spring 2008 – 204 99

  99. How do ants know where to go ? (AS contd.) Ant at point i : ◦ τ ij intensity of pheromone trail from i → j ◦ η ij visibility (constant) from i → j ◦ probability to go to j (simplest version): τ ij η ij p ij = � τ ik η ik k Trail update: ¡ ◦ Old pheromone evaporates partly (0 < ρ < 1) ◦ Ant on route i → j with length ℓ ij spreads q pheromone △ τ ij = q ℓ ij ◦ New pheromone tail τ t +1 = ρ τ t ij + △ τ t ij ij NOS (4311010) – M. Gilli Spring 2008 – 205 AS (contd.) ◦ Search area of the ant corresponds to a discrete set of solutions. ◦ Amount of food is associated with an objective function. ◦ Pheromone trail is modelled with an adaptive memory. 1: Initialize pheromone trail 2: while stopping criteria not met do for all ants do 3: Deposit ant randomly 4: while solution incomplete do 5: Select next element randomly according to pheromone trail 6: end while 7: Compute objective function, Update best solution 8: end for 9: Let a proportion of pheromone evaporate 10: Update pheromone trail (more for better solutions) 11: 12: end while Not so flexible !!! NOS (4311010) – M. Gilli Spring 2008 – 206 Applications: ◦ Travelling salesman problem ◦ Quadratic assignment problem ◦ Job scheduling problem ◦ Graph coloring problem ◦ Sequential ordering References: ◦ Colorni et al. (1992 a ) and Colorni et al. (1992 b ) ◦ Overview on different versions and applications: Bonabeau et al. (1999) ◦ Applications to finance in Maringer (2005) NOS (4311010) – M. Gilli Spring 2008 – note 1 of slide 206 100

Download Presentation
Download Policy: The content available on the website is offered to you 'AS IS' for your personal information and use only. It cannot be commercialized, licensed, or distributed on other websites without prior consent from the author. To download a presentation, simply click this link. If you encounter any difficulties during the download process, it's possible that the publisher has removed the file from their server.

Recommend

NUMERICAL SIMULATION OF SELF-STRUCTURING ANTENNAS BASED ON A GENETIC ALGORITHM OPTIMIZATION SCHEME

NUMERICAL SIMULATION OF SELF-STRUCTURING ANTENNAS BASED ON A GENETIC ALGORITHM OPTIMIZATION SCHEME

NUMERICAL SIMULATION OF SELF-STRUCTURING ANTENNAS BASED ON A GENETIC ALGORITHM OPTIMIZATION SCHEME J.E. Ross * E.J. Rothwell, C.M. Coleman L.L. Nagy John Ross &amp; Associates ECE Dept. Delphi Automotive Systems 350 W 800 N, Suite 317

681 views • 32 slides
Numerical simulation Numerical simulation of dynamics of atmosphere of atmosphere of dynamics

Numerical simulation Numerical simulation of dynamics of atmosphere of atmosphere of dynamics

International conference and young scientists school CITES-2007, Tomsk Lomonosov Moscow State University V. . M M. . Stepanenko Stepanenko, , V D. . N N . . Mikushin Mikushin D Numerical simulation Numerical simulation of dynamics

536 views • 21 slides
Numerical Simulation of Lasers Numerical Simulation of Lasers Christoph P fl aum . p.1/116

Numerical Simulation of Lasers Numerical Simulation of Lasers Christoph P fl aum . p.1/116

Numerical Simulation of Lasers Numerical Simulation of Lasers Christoph P fl aum . p.1/116 Basic Elements of a Laser A laser consists mainly of the following three elements : 1. Laser medium: collection of atoms, molecules, ions or a

1.34k views • 116 slides
to Simulation Optimization to Simulation, Optimization, and back Biology Physics

to Simulation Optimization to Simulation, Optimization, and back Biology Physics

EDA among Sci &amp; Eng fields EDA among Sci &amp; Eng fields CS CS EDA EDA EE EE From Computability to Simulation Optimization to Simulation, Optimization, and back Biology Physics Igor Markov University of Michigan

299 views • 4 slides
Simulation based analysis of human push recovery motions using numerical optimization

Simulation based analysis of human push recovery motions using numerical optimization

, Simulation based analysis of human push recovery motions using numerical optimization malin.schemschat@iwr.uni-heidelberg.de September 1, 2016 malin.schemschat@iwr.uni-heidelberg.de , September 1, 2016 0 / 21 Motivation , adult eldery

533 views • 26 slides
Some Topics in Optimization for Simulation Michael Fu University of Maryland The ACNW

Some Topics in Optimization for Simulation Michael Fu University of Maryland The ACNW

Some Topics in Optimization for Simulation Michael Fu University of Maryland The ACNW OPTIMIZATION TUTORIALS June 13, 2005 Tres Tapas A Soft Overview: Optimization for Simulation (JOC &amp; Annals OR papers) A Little Tutorial:

1.01k views • 54 slides
Systemizing the Solution of Simulation-Driven Optimization Problems Marco Enriquez (joint work

Systemizing the Solution of Simulation-Driven Optimization Problems Marco Enriquez (joint work

Problem Computations Numerical Results Systemizing the Solution of Simulation-Driven Optimization Problems Marco Enriquez (joint work with William Symes) Center for Advanced Aviation Systems Development (CAASD) menriquez@mitre.org November

763 views • 51 slides
Basics of Numerical Optimization: Computing Derivatives Ju Sun Computer Science &amp;

Basics of Numerical Optimization: Computing Derivatives Ju Sun Computer Science &amp;

Basics of Numerical Optimization: Computing Derivatives Ju Sun Computer Science &amp; Engineering University of Minnesota, Twin Cities February 25, 2020 1 / 36 Derivatives for numerical optimization gradient descent Newtons

1.44k views • 105 slides
Tutorial Verification Numerical Simulation Analytic and Numerical Methods in ODEs Further

Tutorial Verification Numerical Simulation Analytic and Numerical Methods in ODEs Further

BENG 221 Lecture 4 Tutorial BENG 221 Mathematical Methods in Bioengineering Overview Example: Ring Oscillator Dynamics Lecture 4 Analytic ODE Solution Numerical Tutorial Verification Numerical Simulation Analytic and Numerical

423 views • 21 slides
AWE Aldermaston UK Scope: 2D and 3D numerical simulation (DNS or LES) of the non-linear growth

AWE Aldermaston UK Scope: 2D and 3D numerical simulation (DNS or LES) of the non-linear growth

REVIEW OF NUMERICAL SIMULATION OF MIXING DUE TO RAYLEIGH-TAYLOR AND RICHTMYER-MESHKOV INSTABILITIES David Youngs AWE Aldermaston UK Scope: 2D and 3D numerical simulation (DNS or LES) of the non-linear growth of Rayleigh-Taylor and

740 views • 61 slides
Outline Narcisse Ngada DESY, MKK 1) What is simulation ? 14.05.2014 2) Why simulation ? 3)

Outline Narcisse Ngada DESY, MKK 1) What is simulation ? 14.05.2014 2) Why simulation ? 3)

Simulation Outline Narcisse Ngada DESY, MKK 1) What is simulation ? 14.05.2014 2) Why simulation ? 3) Principles of simulation 4) Types of simulation Analog simulation Numerical simulation 5) Conclusion 1. What is simulation ? 2.

1.08k views • 54 slides
Examples of Numerical Modeling and Simulation for Automotive EMC and RF Applications Marco

Examples of Numerical Modeling and Simulation for Automotive EMC and RF Applications Marco

Examples of Numerical Modeling and Simulation for Automotive EMC and RF Applications Marco KLINGLER &amp; Salah BENHASSINE PSA Peugeot Citron Confrence La simulation en CEM et Hyperfrquences , Salon Microwave &amp; RF, Paris, 3-5

586 views • 40 slides
Optimization and Simulation Introduction to simulation Michel Bierlaire Transport and Mobility

Optimization and Simulation Introduction to simulation Michel Bierlaire Transport and Mobility

Optimization and Simulation Introduction to simulation Michel Bierlaire Transport and Mobility Laboratory School of Architecture, Civil and Environmental Engineering Ecole Polytechnique F ed erale de Lausanne M. Bierlaire (TRANSP-OR ENAC

864 views • 32 slides
Glazunov A.V. Institute of Numerical Mathematics RAS glazunov@inm.ras.ru Large-eddy simulation

Glazunov A.V. Institute of Numerical Mathematics RAS glazunov@inm.ras.ru Large-eddy simulation

Glazunov A.V. Institute of Numerical Mathematics RAS glazunov@inm.ras.ru Large-eddy simulation of shear-driven turbulence using dynamic SFS closure and high-order numerical scheme Synoptical variations Boundary-Layer flows DNS requirement

715 views • 44 slides
Numerical Simulation of the Stress-Strain Behavior of Ni-Mn-Ga Shape Memory Alloys Marcel Arndt

Numerical Simulation of the Stress-Strain Behavior of Ni-Mn-Ga Shape Memory Alloys Marcel Arndt

Institute for Numerical Simulation, University of Bonn Numerical Simulation of the Stress-Strain Behavior of Ni-Mn-Ga Shape Memory Alloys Marcel Arndt arndt@ins.uni-bonn.de Institute for Numerical Simulation Rheinische

677 views • 22 slides
Optimization and Simulation Constrained optimization Michel Bierlaire Transport and Mobility

Optimization and Simulation Constrained optimization Michel Bierlaire Transport and Mobility

Optimization and Simulation Constrained optimization Michel Bierlaire Transport and Mobility Laboratory School of Architecture, Civil and Environmental Engineering Ecole Polytechnique F ed erale de Lausanne M. Bierlaire (TRANSP-OR ENAC

1.02k views • 72 slides
A Topology-Aware Performance Monitoring Tool for Shared Resource Management in Multicore Systems

A Topology-Aware Performance Monitoring Tool for Shared Resource Management in Multicore Systems

A Topology-Aware Performance Monitoring Tool for Shared Resource Management in Multicore Systems TADaaM Team - Nicolas Denoyelle - Brice Goglin - Emmanuel Jeannot August 24, 2015 1. Context/Motivations 2. Fast presentation of the tool 3.

446 views • 29 slides
SAPM Overview In this lecture we review the topics we have covered this Semester Summary

SAPM Overview In this lecture we review the topics we have covered this Semester Summary

SAPM Overview In this lecture we review the topics we have covered this Semester Summary semester, focusing on what I consider the most important points to remember. Dr. James A. Bednar The lecture slides on each topic, coupled with the

71 views • 6 slides
Practical Bioinformatics Mark Voorhies 4/28/2017 Mark Voorhies Practical Bioinformatics

Practical Bioinformatics Mark Voorhies 4/28/2017 Mark Voorhies Practical Bioinformatics

Practical Bioinformatics Mark Voorhies 4/28/2017 Mark Voorhies Practical Bioinformatics Pearson distances Pearson similarity N i ( x i x offset )( y i y offset ) s ( x , y ) = N N i ( x i x offset ) 2 i ( y i

502 views • 28 slides
Lists and Tuples Ali Taheri Sharif University of Technology Spring 2019 Outline 1. List Basics

Lists and Tuples Ali Taheri Sharif University of Technology Spring 2019 Outline 1. List Basics

Fu Fundamentals of Pr Programming (Py Python) Lists and Tuples Ali Taheri Sharif University of Technology Spring 2019 Outline 1. List Basics 2. List Operations 3. Lists are Mutable 4. Iterating over Lists 5. List Methods 6. Matrices

375 views • 16 slides
Registres ANOCOR (ANOmalies CORonaires congnitales) Xavier HALNA du FRETAY Unit

Registres ANOCOR (ANOmalies CORonaires congnitales) Xavier HALNA du FRETAY Unit

Registres ANOCOR (ANOmalies CORonaires congnitales) Xavier HALNA du FRETAY Unit Cardiologique de la Reine Blanche, SARAN D CLARATION DE LIENS DINTRT AVEC LA PRSENTATION Intervenant : Halna du Fretay Xavier, SARAN, FRANCE Je

353 views • 23 slides
robuBOX-Kompa : Open Source software for companion robots May 18 nd , 2010 Franois HIRIGOYEN

robuBOX-Kompa : Open Source software for companion robots May 18 nd , 2010 Franois HIRIGOYEN

robuBOX-Kompa : Open Source software for companion robots May 18 nd , 2010 Franois HIRIGOYEN Consuelo GRANATA Summary Summary 1. About us 2. Kompa - a closer look at the platform 3. What is PURE? 4. How is it implemented? 5. A

800 views • 45 slides
Disclosures Co-author patent application regarding Sinusitis diagnostics and treatments

Disclosures Co-author patent application regarding Sinusitis diagnostics and treatments

Surgery of the Frontal Sinus Steven D. Pletcher MD Assistant Professor Department of Otolaryngology Head and Neck Surgery University of California, San Francisco Disclosures Co-author patent application regarding Sinusitis

954 views • 28 slides
PROMYS FOR T EACHERS (P F T) Ryota Matsuura St. Olaf College March 24, 2011 A RNOLD R OSS ONCE

PROMYS FOR T EACHERS (P F T) Ryota Matsuura St. Olaf College March 24, 2011 A RNOLD R OSS ONCE

P ANEL : S ECONDARY T EACHER P ROGRAMS Gail Burrill Al Cuoco Ryota Matsuura March 24, 2011 P ANEL OVERVIEW We will describe three distinct but closely related mathematical experiences for secondary teachers: PROMYS for Teachers (Ryota) Focus

800 views • 39 slides

More recommend