# Optimization for Stochastic Simulations

This questions stems from my original question here. I have a system, where each 'run' is a simulation of agents switching between tasks (or states) as in a Markov Chain. I would like to find ways to optimize the Transition Probability Matrix (TPM) to achieve a desired distribution of the agents among the states (which is determined by each run of the simulation) in the steady state.

What exactly do I call such a system, and where can I learn more about working with these systems (names of topics, keywords, authors etc.)?

I am aware of the analytical methods and relations for arriving at the exact values in the TPM.

• Implemented this using robustified damped BFGS Quasi-Newton method w/ finite difference gradients. Wanted to use Reverse Mode Automatic Differentiation of objective function, but tool I was using, ADiMat, didn't work correctly on this (used it successfully for other stochastic simulations). Automatic Differentiation interchanges differentiation and expectation, which works when objective is continuous per replication. I've even solved it adding a constraint that Transition Matrix is symmetric positive semidefinite, which arises in some problems arxiv.org/pdf/1301.4055.pdf It works!! Oct 30, 2022 at 14:09
• I did not use evolutionary algorithm or any algorithm named for purported animal behavior. I used derivative based methods, robustified to deal with the noise, and to determine when to adjust the number of replications performed, and to assess when the solition has either converged to an optimum, or is as good as it gets if the max allowed number of replications per iteration is already being used - all of this based on stochasticized KKT optimality criteria. Note:Automatic Differentiation, allows you to use stochastic simulation to evaluate the gradient,without using finite differences. Oct 30, 2022 at 14:18
• If you have a specific problem instance you are interested in, feel free to edit your question to include it. You can also include additional constraints, linear or nonlinear, including Linear or Nonlinear SDP (matrix constrained to be positive semidefnite) constraints; and my program can handle it. I have implemented a last squares objective (which minimizes sum of squares of difference between stationary distribution corresponding the optimized Transition Matrix, and the target stationary distribution), but can change that to suit your needs. Oct 30, 2022 at 14:28
• Tying other variants: Assume adversarial game. Someone trying to maximize, not minimize, squared difference to target stationary distribution., Objective function very far from being convex, anywhere, so BFGS Quasi-Newton likely to perform poorly, even without noise. So I use noise resistant variant of SR1, which allows eigenvalues of Hessian to be any mixture of positive, negative, zero. Transition Matrix also constrained to be symmetric positive semidefinite. Due to non-convex objective, this is a non-convex Nonlinear SDP, a formidable beast, even without noise,. But I'm giving it a go. Oct 30, 2022 at 17:09
• And it solved the non-convex Nonlinear SDP with stochastic simulation objective function like a champ!! No prob, Bob. Due to most variables quickly moving to a bound, which turned out to be optimal, it was easier to solve than the Linear SDP problem. Nov 3, 2022 at 2:57

Simulation Optimization is the name of the field you are looking for. It is closely related to, and overlaps with. Machine Learning, such as training Deep Learning models.

In Simulation Optimization, there is no such thing as training on all the data, because there are an infinite number of data points which could be generated, not only at different values of the optimization variables, but the number of replications increasing to infinity.

There are a number of general approaches, including stochastic gradient descent (originally called Kiefer-Wolfowitz when it was invented 70 years ago, although that was the finite difference version, and only in one dimension), Gradient based methods, to include Quasi-Newton, or even Newton (simulating for Hessian).

Any algorithm which uses averages of sample for objective function and gradients (or Hessian) can be called Sample Average Approximation (SAA). The term is generic enough that it can refer to anything. But SAA is often used as a synonym for any algorithm using common random numbers, which is also called Sample Path Optimization, although that too can mean almost anything.

It is noteworthy, especially in conjunction with derivative based methods, but also for zeroth order (Objective function only) methods, use of common random numbers can be extremely important for performance. If calculating a finite difference derivative at two nearby points, if the errors are independent, the variances add, and get divided by a small number, which creates a huge variance in the finite difference derivative; that's not a good thing. However, if common random numbers are used, and particularity if there is good synchronization in the implementation so that very high positive correlation is achieved at nearby points, the individual data point errors might cancel out to a large extent, and finite difference derivatives could potentially be more accurate than evaluation of the objective at an individual point Similarly for differencing gradients from adjacent iterations in order to feed a Quasi-Newton Hessian update - this concept seems foreign to many Machine Learning people.

Surrogate models (sometime known as Response Surface Methodology) fit a model to objective values at sampled points, and then optimize the surrogate function; then additional points ere sampled to improve the model, then re-optimized, etc. One manifestation is fitting a Gaussian Process, then optimizing that. (a form of Bayesian Optimization) Stochastic Kriging is another term for essentially the same thing. The uncertainty in the fitted model expands between sampled data points. There are also a variety of smoothing methods, involving some form of splines, or otherwise.

Evolutionary algorithms and various "animal" algorithms are used by some people, but I don't have a very positive view of them. Perhaps they have their use. You seem to have a handle on those, so I won't further elaborate.

Going back to derivative based methods, there are several possibilities. Finite difference gradient are always applicable if the objective function is theoretically differentiable. A method called Likelihood Ratio Method (it was my question over 40 years ago to the then student who developed LRM, of how you could simulate for a derivative, which I wanted for use in optimization, that stimulated him to invent the method) is fairly broadly applicable to directly simulate for derivatives, but is usually not very fast (has high variance), and not used that widely - it basically involves differentiating under the expectation, but instead of differentiating the (unintegrated) objective function, i.e., the integrand, it differentiates the probability measure driving the simulation) A closely related method, known as week derivatives, is essentially a more modern version of Likelihood Ratio Method, but as of yet, no one has created an Automatic Differentiation (AD) like machinery to make it "easy" to use on complicated models - I think that it might be doable by marrying the expertise of an AD expert to modify existing AD tools to incorporate different end node rules needed for weak derivatives (there might be some theoretical applicability issues though).

If the objective function is continuous (smooth) per sample path, not just in expected value, then it is generally valid to interchange differentiation and expectation (simulation evaluates expectation, even when estimating probabilities). This is known as Infinitesimal Perturbation Analysis (IPA). The objective function is differentiated, hen simulated. Automatic Differentiation (a.k.a. Algorithmic Differentiation) can be used to directly simulate for the gradient, especially using the so-called Reverse Mode of Automatic Differentiation, which computes the entire gradient, no matter how many components, at a fixed overhead ratio to just the objective function (presuming three is enough memento). The Reverse Mode of AD applied to Neural Networks is better known as Backpropagation. Sometimes Backpropagation is used to refer just to the Reveres Mode AD, while sometime it is used sued to refer to the whole stochastic gradient descent algorithm using Reverse Mode AD.

IPA using the Reverse Mode of Automatic Differentiation can compute (simulate) an entire gradient "at once". By contrast, LRM calculates the gradient one component at a time, by separately simulating with respect to the probability measure (partially) differentiated with respect to one component at a time. This can be a big disadvantage relative to IPA in models having a large number of optimization variables.

Complex Step Differentiation is a method of doing finite differences using complex numbers. Despite looking very different from Automatic Differentiation, it essentially amounts to the same thing (once the simulation is suitably modified to handle complex inputs), and achieves the same accuracy as Automatic Differentiation,. Unfortunately, its scope of applicability is the same as IPA, not that of finite differences, i.e., it requires each sample path be differentiable, not just in expectation. Further, Complex Step Differentiation really corresponds to the Forward Mode of Automatic Differentiation; and therefore, computes only one component of the gradient at a time. This places it at a disadvantage to RM Automatic Differentiation when the number of optimization variables sis large, but might be advantageous for a small number of variables, due to avoiding the overhead of RM Automatic Differentiation.

If some or all of the optimization variables are discrete, that's a whole different ball game which I won't get into here.

• Props for a very thorough answer.
– prubin
Nov 1, 2022 at 18:40
• Weak Derivatives were developed after LRM, and kind of LRM on steroids. MIchael C. Fu "What you should know about simulation and derivatives" onlinelibrary.wiley.com/doi/abs/10.1002/nav.20313 is great paper, especially table on p. 729 listing derivatives of probability distributions for IPA, LRM and WD (weak derivatives). Main merit to LRM or WD is when the sample paths is not continuous per realization (sample path), and can't easily be made so. For large number of variables, as in Deep Learning, IPA, implemented by Reverse Mode AD is the hands down winner, hence its popularity. Nov 3, 2022 at 14:24
• As I suggested, implementing Weak Derivatives via modified Automatic Differentiation (AD) (at least Forward Mode, not sure if Reverse Mode is doable) would be great contribution - perhaps nice Ph.D. project. You'd have to work with an expert in AD and perhaps modify existing AD SW to insert derivatives of probability measures at modes instead of 'regular" derivatives. That could make WD viable for widespread use There might be some regularity conditions to worry about, which you can try to address, or just differentiate those bad boys, regularity conditions be damned. and let god sort it out Nov 3, 2022 at 14:33
• In early 1982, at the same time I was inquiring with Peter Glynn, who developed LRM, and unbeknownst to me, Ho, Cao, Cassandras were developing IPA, I independently "invented" 'IPA - I didn't call it that, just called it interchanging differentiation and expectation. But I dismissed it as a practical solution for optimization because I thought the technical conditions required for interchange wouldn't be satisfied very often in real-world simulations, and even if it were, would be almost impossible to verify in complicated models. It turns out I was wrong - it is very widely applicable. Nov 3, 2022 at 23:50
• Before 1982, no one thought about simulating for derivatives. I recognized & wanted to fill that gap so I could use derivative-based methods for simulation optimization, which no one contemplated either, back then.. PWG told me a week after I inquired that he figured out how to simulate for derivatives. All he said was "Likelihood Ratio". I asked him how that was used to simulate for derivatives. He didn't say, he just smirked and said "Likelihood Ratio". It wasn't until 4 years later when he published the first paper on LRM that I finally found out what Likelihood Ratio meant in that context Nov 4, 2022 at 0:06