Two Topics in Parametric Integration Applied to Stochastic - - PowerPoint PPT Presentation

Two Topics in Parametric Integration Applied to Stochastic Simulation in Industrial Engineering

Jeremy Staum

Department of Industrial Engineering and Management Sciences Northwestern University

September 15th, 2014

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

Outline

1 Simulation Metamodeling

introduction and overview

2 Multi-Level Monte Carlo Metamodeling

with Imry Rosenbaum http: //users.iems.northwestern.edu/~staum/MLMCM.pdf

3 Generalized Integrated Brownian Fields for Simulation

Metamodeling

with Peter Salemi and Barry L. Nelson http://users.iems.northwestern.edu/~staumGIBF.pdf

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

MCQMC / IBC Application Domain

Industrial Engineering & Operations Research using math to analyze systems and improve decisions Stochastic Simulation: production, logistics, financial, . . . integration: µ = E[Y ] =

• Y (ω) dω

parametric integration: approx. µ def. by µ(x) = E[Y (x)]

• ptimization: min{µ(x) : x ∈ X}

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

What is Stochastic Simulation Metamodeling?

Stochastic simulation model example Fuel injector production line System performance measure µ(x) = E[Y (x)]

x: design of production line Y (x): number of fuel injectors produced

Simulating each scenario (20 replications) takes 8 hours Stochastic simulation metamodeling Simulation output ¯ Y (xi) at xi, i = 1, . . . , n Predict µ(x) by ˆ µ(x), even without simulating at x ˆ µ(x) is usually a weighted average of ¯ Y (x1), . . . , ¯ Y (xn)

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

Overview of Multi-Level Monte Carlo (MLMC)

Error in Stochastic Simulation Metamodeling prediction ˆ µ(x) =

k

• i=1

wi(x) ¯ Y (xi) variance: Var[ˆ µ(x)] caused by variance of simulation output interpolation error (bias): E[ˆ µ(x)] = µ(x)

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

Main Idea of Multi-Level Monte Carlo

Ordinary Monte Carlo to reduce variance: large number n of replications per simulation run (design point) to reduce bias: large number k of design points (fine grid) very large computational effort kn Multi-Level Monte Carlo to reduce variance: coarser grids, many replications each to reduce bias: finer grids, few replications each less computational effort / better convergence rate

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

Our Contributions

Theoretical mix and match or expand (from Heinrich’s papers): derive desired conclusions under desired assumptions to suit IE goals and applications Practical algorithm design (based on Giles) Experimental show how much MLMC speeds up realistic examples in IE

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

Heinrich (2001): MLMC for Parametric Integration

Assumptions Approximate µ given by µ(x) =

• Ω Y (x, ω) dω over x ∈ X.

X ⊆ Rd and Ω ⊆ Rd2: bounded, open, Lipschitz boundary. With respect to x, Y has weak derivatives up to rth order. Y and weak derivatives are Lq-integrable in (x, ω). Sobolev embedding condition: r/d > 1/q. Measure error as (

• Ω ˆ

µ − µp

q dω)1/p, where p = min{2, q}.

Conclusion: There is a MLMC method with optimal rate. MLMC attains the best rate of convergence in C, the number of evaluations of Y . The error bound is proportional to C −r/d if r/d < 1 − 1/p C 1/p−1 log C if r/d = 1 − 1/p C 1/p−1 if r/d > 1 − 1/p.

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

Assumptions

Smoothness: assume r = 1 Stock option, Y (x, ω) = max{xR(ω) − K, 0} Queueing: waiting time Wn+1 = max{Wn + Bn − An+1, 0} Inventory: Sn = min{In + Pn, Dn}, In+1 = In − Sn Parameter Domain Assume X ⊂ Rd is compact (not open).

Heinrich and Sindambiwe (1999), Daun and Heinrich (2014)

If X were open, we would have to extrapolate. No need to approximate unbounded µ near a boundary of X. Domain of Integration Ω ⊆ Rd2 is not important; d2 does not appear in theorem.

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

Changing Perspective

Measure of Error Use p = q = 2 to get Root Mean Integrated Squared Error

• Ω
• X (ˆ

µ(x) − µ(x))2 dx dω 1/2

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

Changing Perspective

Measure of Error Use p = q = 2 to get Root Mean Integrated Squared Error

• Ω
• X (ˆ

µ(x) − µ(x))2 dx dω 1/2 Sobolev Embedding Criterion with r = 1, q = 2 r/d > 1/q becomes 1/d > 1/2, i.e. d = 1!??

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

Changing Perspective

Measure of Error Use p = q = 2 to get Root Mean Integrated Squared Error

• Ω
• X (ˆ

µ(x) − µ(x))2 dx dω 1/2 Sobolev Embedding Criterion with r = 1, q = 2 r/d > 1/q becomes 1/d > 1/2, i.e. d = 1!?? Why We Don’t Need the Sobolev Embedding Condition Assume the domain X is compact. Assume Y (·, ω) is (almost surely) Lipschitz continuous. Conclude Y (·, ω) is (almost surely) bounded.

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

Our Assumptions

On the Stochastic Simulation Metamodeling Problem X ⊂ Rd is compact Y (x) has finite variance for all x ∈ X |Y (x′, ω) − Y (x, ω)| ≤ κ(ω)x′ − x, a.s., and E[κ2] < ∞. On the Approximation Method and MLMC Design ˆ µ(x) = N

i=1 wi(x) ¯

Y (xi) where each wi(x) ≥ 0 and Total weight on points xi far from x gets close to 0. Total weight on points xi near x gets close to 1. Thresholds for “far”/“‘near” and “close to” are O(N−1/2φ) as number N of points increases. Examples: piecewise linear interpolation on a grid; nearest-neighbors, Shepard’s method, kernel smoothing

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

Approximation Method Used in Examples

Kernel Smoothing ˆ µ(x) =

N

• i=1

wi(x) ¯ Y (xi) weight wi(x) is

0 if xi is outside the cell containing x

• therwise, proportional to

exp(−x − xi)

weights are normalized to sum to 1

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

Our Conclusions

MLMC Performance As number N of points used in a level increases, Errors due to bias and refinement variance are like O(N−1/φ). Example: nearest-neighbor approximation on grid, φ = d/2 Computational Complexity (based on Giles 2013) To attain RMISE < ǫ, the required number of evaluations of Y is O(ǫ−2(1+φ)) for standard Monte Carlo and for MLMC it is O(ǫ−2φ) if φ > 1 O((ǫ−1(log ǫ−1))2) if φ = 1 O(ǫ−2) if φ < 1.

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

Sketch of Algorithm (based on Giles 2008)

Goal: add levels until target RMISE < ǫ is achieved.

1 INITIALIZE level ℓ = 0. 2 SIMULATE at level ℓ: 1

Run level ℓ simulation experiment with M0 replications.

2

Observe sample variance of simulation output.

3

Choose number of replications Mℓ to control variance; run more replications if needed.

3 TEST CONVERGENCE: 1

Use Monte Carlo to estimate the size of the refinement ∆ˆ µℓ,

• X (∆ˆ

µℓ(x))2 dx.

2

If refinements are too large compared to target RMISE, increment ℓ and return to step 2.

4 CLEAN UP: Finalize number of replications M0, . . . , Mℓ to

control variance; run more replications at each level if needed.

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

Asian Option Example, d = 3

MLMC up to 150 times better than standard Monte Carlo

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

Inventory System Example, d = 2

MLMC was 130-8900 times better than standard Monte Carlo

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

Conclusion on Multi-Level Monte Carlo

Celebration Multi-Level Monte Carlo works for typical IE stochastic simulation metamodeling too! Future Research Handle discontinuities in simulation output. Combine with good experiment designs. Grids are not good in high dimension.

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

Introduction: Generalized Integrated Brownian Field

Kriging / Interpolating Splines Pretend µ is a realization of a Gaussian random field M with mean function m and covariance function σ2. Kriging predictor: ˆ µ(x) = m(x) + σ2(x)Σ−1(¯ Y − m) = m(x) +

• i

βiσ2(x, xi)

σ2(x) is a vector with ith element σ2(x, xi) Σ is a matrix with i, jth element σ2(xi, xj) ¯ Y − m is a vector with ith element ¯ Y (xi) − m(xi)

Stochastic Kriging / Smoothing Splines ˆ µ(x) = m(x) + σ2(x)(Σ + ❈)−1(¯ Y − m) = m(x) +

• i

˜ βiσ2(x, xi) ❈ = covariance matrix of noise, estimated from replications

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

Radial Basis Functions vs. Integrated Brownian Field

Radial Basis Functions Basis Functions from r-Fold Integrated Brownian Field

(d) r = 0 (e) r = 1 (f) r = 2

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

Response Surfaces in IE Stochastic Simulation

(g) Credit Risk (h) Inventory

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

r-Integrated Brownian Field Br

Covariance function / reproducing kernel σ2(x, y) =

d

• i=1

1 (r!)2 1 (xi − ui)r

+(yi − ui)r + dui

Inner product f , g =

• (0,1)d(f ([r···r])(u))(g([r···r])(u)) du

Space Tensor product of Sobolev Hilbert space Hr(0, 1) with boundary conditions f (j)(0) = 0 for j = 0, . . . , r What’s missing? polynomials of degree r

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

Removing Boundary Conditions: d = 1

Generalized integrated Brownian motion Xr(x) =

r

• k=0
• θkZk

xk k! +

• θr+1Br(x)

Covariance function / reproducing kernel σ2(x, y) =

r

• k=0

θk xkyk (k!)2 + θr+1 1 (x − u)r

+(y − u)r +

(r!)2 du Sobolev space Hr(0, 1) , no boundary conditions Inner product f , g =

r

• k=0

1 θk (f (k)(u))(g(k)(u)) + 1 θr+1 1 (f (r)(u))(g(r)(u)) du

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

Multidimensional, Without Boundary Conditions

Tensor-Product RKHS with Weights Example of reproducing kernel for d = 2, r = 1 K(x, y) = θ00 + θ10x1y1 + θ20(x1 ∧ y1) + θ01x2y2 + θ02(x2 ∧ y2) +θ11x1x2y1y2 + θ12x1y1(x2 ∧ y2) +θ21(x1 ∧ y1)x2y2 + θ22(x1 ∧ y1)(x2 ∧ y2) In general, one weight for each of d

i=1(ri + 2) subspaces.

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

Multidimensional, Without Boundary Conditions

Tensor-Product RKHS with Weights Example of reproducing kernel for d = 2, r = 1 K(x, y) = θ00 + θ10x1y1 + θ20(x1 ∧ y1) + θ01x2y2 + θ02(x2 ∧ y2) +θ11x1x2y1y2 + θ12x1y1(x2 ∧ y2) +θ21(x1 ∧ y1)x2y2 + θ22(x1 ∧ y1)(x2 ∧ y2) In general, one weight for each of d

i=1(ri + 2) subspaces.

Generalized Integrated Brownian Field Covariance function / reproducing kernel σ2(x, y) =

d

• i=1

ri

• k=0

θi,k xk

i yk i

(k!)2 + θi,ri+1 1 (xi − ui)ri

+(yi − ui)ri +

(ri!)2 dui

• In general, number of weights is d

i=1(ri + 2).

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

Our Contributions

more parsimonious parametrization

makes maximum likelihood estimation easier and MLE search for r1, . . . , rd

GIBF has Markov property

d = 1: proof d > 1: conjecture

IE simulation examples

stochastic and deterministic simulation standard and nonstandard information

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

Credit Risk Example, d = 2

Experiment design: 63 Sobol’ points, predictions in a smaller square Factor by which MISE decreased using (1,1)-GIBF Number of replications ∞ 100 25 Noise level none low medium Without gradient estimates 94 111 120 With gradient estimates 81 83 69

(i) Credit risk surface (j) Gaussian (k) (1, 1)-GIBF

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

Conclusion on Generalized Integrated Brownian Field

Emancipating Simulation Metamodeling from Geostatistics a new covariance function for kriging, designed for simulation metamodeling in engineering Superior Practical Performance 4-120 times better than Gaussian covariance function in 2-6 dimensional examples with or without gradient information

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering

Thank You!

Jeremy Staum Topics in Simulation Metamodeling in Industrial Engineering