How can simulations be speeded up by orders of magnitude

Monte Carlo simulation are often used as the performance evaluation tool of choice for the analysis and optimization of complex human-made systems such as manufacturing plants, communication networks, logistic-command-control systems and other civilian/military systems. Mathematically, we are interested in the evaluation and optimization of the performance of the system represented by Min θ ∈ Θ J(θ)-E[L(θ,ξ)] (1) where Θ, the search space, is an arbitrary, huge, but finite set; θ the design alternatives; J, the performance criterion which is the expectation of L, the sample performance, as a functional of θ and ξ, the randomness in the systems. This definition is sufficiently general to include all sorts of performance analysis and synthesis problems of interest. There are two aspects to this problem: the stochastic part and the optimization part. The stochastic aspect has to do with the problem of performing numerical expectation since the functional L(x(t;θ,ξ)) is available only in the form of a complex calculation or simulation. The standard approach is to estimate E[L(θ,ξ)] by Monte Carlo sampling, i.e., N E[L(θ,ξ)]≃1/NΣ N i=1 L(θ,ξ i )]?J(θ) est (2) where ξ i represents the ith sample of the randomness. However, the ultimate accuracy (confidence interval) of such an estimate cannot be improved upon faster than 1/(N) 1/2 , the result of averaging i.i.d. noise. Here each sample of L(θ,ξ i ) requires one simulation run which for complex problems can take seconds or minutes even on the fastest computer not to mention the setup cost of programming the simulation. Millions or billions of required samples of L(θ,ξ i ) for different value of θ may become an infeasible burden. Yet there is no getting away from this fundamental limit of 1/(N) 1/2 , i.e. one order of magnitude increase in accuracy requires two orders of magnitude increase in sampling cost. The second limitation has to do with the optimization part. Traditionally, analysis again plays an important part in numerical optimization whether it is convergence, gradient computation, or general hill climbing. Again tools, such as perturbation analysis (PA), take advantage of the structure and real variable nature of the problem to devise effective algorithms for optimization. The best known example is of course linear programming. However, when we exhaust such advantages and the remaining problem becomes structureless and the search space Θ is totally arbitrary, we encounter the limitation introduced by combinatorial explosion which we shall call the NP-hard limitation'. In such case, more or less blind search is the only alternative. While it is often easy to prove probabilistic convergence of random search, it is very hard to improve the search efficiency. An example will illustrate this: Suppose Θ is of size 10 10 which is small by combinatorial standards. We are allowed to take 10,000 samples of J(θ) uniformly in Θ. What is the probability that at least one of these sample will belong in the top-50, top-500, or top-5000 of the designs of Θ? Simple calculation yields, Prob(at least one of 10,000 sample is in top-k)=1-(1-k/10 10 ) 10,000 (3) which turns out to be 0.00004999, 0.0004998, and 0.004988 respectively, i.e., forget it! This means that for simulation based optimization problem of any complexity, simple random search is not an effective approach - too much computation for too little return.