The Simulation Optimization Problem

General continuous SO problem:

\[\min_{\mathbf{x} \in X \subseteq \mathbb{R}^d} \; E[G(\mathbf{x})]\]

where:

• \(G : \mathbf{x} \mapsto \mathbb{R}\) is a black-box objective, observable only through a stochastic simulation oracle.
• At any \(\mathbf{x}\), only a noisy estimate \(\hat{G}(\mathbf{x})\) is available.
• The feasible region \(X\) may be continuous, integer-ordered, or mixed.

Practical examples: system configuration parameters, staffing levels, inventory policies, service rates, resource allocation vectors.

Where CE Fits in the SO Landscape

CE targets global search on continuous (or mixed) domains.

Key properties that make CE attractive for simulation optimization:

• Derivative-free — uses only function evaluations, no gradient information needed.
• Noise-tolerant — the population-based update naturally averages out stochastic noise.
• Self-tuning — the sampling distribution automatically narrows around good regions.
• Flexible — the parametric family \(\{f(\cdot;\mathbf{v})\}\) can be chosen to match the problem structure.

CE Origins and Scope

The CE method was originally developed for rare event probability estimation (Rubinstein, 1997).

It was quickly recognized that the same adaptive sampling framework could be redirected toward combinatorial and continuous optimization (Rubinstein, 1999).

In this presentation:

• Rare event simulation is mentioned only as motivation for the underlying mathematics.
• The focus is entirely on CE as a simulation optimization algorithm.
• All KSL implementation details refer to the optimization use case.

The Associated Stochastic Problem (ASP)

Key idea: transform the optimization problem into a probability estimation problem.

For the maximization problem \(\gamma^* = \min_{\mathbf{x} \in X} G(\mathbf{x})\), define for any threshold \(\gamma\):

\[\ell(\gamma) = \mathbb{P}_{\mathbf{u}}(G(\mathbf{X}) \leq \gamma) = \mathbb{E}_{\mathbf{u}} \, I_{\{G(\mathbf{X}) \leq \gamma\}}\]

where \(\mathbf{X}\) is drawn from a parametric family of densities \(f(\cdot; \mathbf{u})\) on \(X\).

The connection to optimization:

• When \(\gamma = \gamma^*\) and \(f(\cdot;\mathbf{u})\) is uniform, \(\ell(\gamma^*)\) is very small — a “rare event.”
• The CE method finds a reference parameter \(\mathbf{v}^*\) such that \(f(\cdot;\mathbf{v}^*)\) concentrates mass near \(\mathbf{x}^*\).
• As \(\mathbf{v}_t \to \mathbf{v}^*\), the sampling distribution degenerates to a point mass at the optimum.

The Parametric Family and Reference Parameter

Choose a parametric family \(\{f(\cdot;\mathbf{v}), \, \mathbf{v} \in V\}\) on \(X\).

• For continuous problems: product of independent Normal distributions, one per dimension.
• For discrete problems: product of independent Bernoulli distributions.
• The parameter vector \(\mathbf{v}\) (called the reference parameter) fully describes the current sampling distribution.

Goal of the CE algorithm: construct a sequence \(\hat{\mathbf{v}}_0, \hat{\mathbf{v}}_1, \hat{\mathbf{v}}_2, \ldots\) such that:

1. The threshold \(\hat{\gamma}_t\) rises toward \(\gamma^*\) at each iteration.
2. \(f(\cdot;\hat{\mathbf{v}}_t)\) concentrates progressively more mass near \(\mathbf{x}^*\).
3. The sequence terminates when the distribution has collapsed — i.e., it has converged to a near-degenerate distribution.

The Two-Phase Loop

Every CE iteration consists of exactly two phases:

repeat until convergence:

  Phase 1 — GENERATE
    Draw N candidates X_1, ..., X_N ~ f(·; v_{t-1})
    Evaluate G(X_i) for each i via the simulation oracle

  Phase 2 — UPDATE
    Identify the elite set: top ρ·N performers
    Re-fit v_t by maximizing the log-likelihood over the elite set
    (Equivalent to computing the MLE of v restricted to the elites)

This is the entire algorithm. All complexity lies in:

• choosing the family \(f(\cdot;\mathbf{v})\),
• deriving the analytic update formula, and
• specifying when to stop.

The Kullback-Leibler Divergence — One Slide

The update formula is derived by minimizing the Kullback-Leibler (KL) divergence between the ideal sampling distribution \(g^*\) and \(f(\cdot;\mathbf{v})\):

\[D(g^*, f(\cdot;\mathbf{v})) = \int g^*(\mathbf{x}) \ln \frac{g^*(\mathbf{x})}{f(\mathbf{x};\mathbf{v})} \, d\mathbf{x}\]

Minimizing \(D\) with respect to \(\mathbf{v}\) is equivalent to solving:

\[\max_{\mathbf{v}} \; \mathbb{E}_{\mathbf{u}} \left[ I_{\{G(\mathbf{X}) \leq \gamma\}} \ln f(\mathbf{X}; \mathbf{v}) \right]\]

Practical consequence: the update step is a weighted maximum likelihood estimation where:

• Only the elite samples (those with \(G(\mathbf{X}_i) \leq \hat{\gamma}_t\)) contribute.
• The update formula is analytic for distributions in the natural exponential family (Normal, Bernoulli, Poisson, etc.).

Algorithm 3.3 — Pseudocode

Initialize: v_0 = u (e.g., problem midpoints for means), t = 1

repeat:

  Step 1.  Generate X_1, ..., X_N  ~  f(·; v_{t-1})

  Step 2.  Compute performances S(X_i) for i = 1, ..., N
           Sort: S_(1) ≤ S_(2) ≤ ... ≤ S_(N)
           Set elite threshold:  γ_t = S_(⌈(ρ)N⌉)   [the (ρ)-quantile]

  Step 3.  Elite set:  E_t = { X_i : S(X_i) ≥ γ_t }

  Step 4.  Update:  v_t = argmax_v  Σ_{i ∈ E_t} ln f(X_i ; v)
                         (analytic MLE restricted to elite samples)

  Step 5.  Apply smoothing:  v_t ← α v_t + (1-α) v_{t-1}

  Step 6.  Check stopping criterion; if met, stop

           t ← t + 1

Output: \(\hat{\mathbf{v}}_T\) — the parameters represent the recommended solution \(\mathbf{x}^* \approx \boldsymbol{\mu}_T\).

The Adaptive Threshold \(\hat{\gamma}_t\)

At each iteration, \(\hat{\gamma}_t\) is the empirical \((\rho)\)-quantile of the sample performances:

\[\hat{\gamma}_t = S_{(\lceil (\rho)N \rceil)}\]

What this achieves:

• \(\hat{\gamma}_t\) rises from iteration to iteration as the sampling distribution improves.
• The elite set always contains approximately \(\lfloor \rho N \rfloor\) samples — enough to re-fit \(\mathbf{v}_t\) reliably.
• The sequence \(\hat{\gamma}_1 \leq \hat{\gamma}_2 \leq \cdots \leq \hat{\gamma}_T\) tracks progress toward \(\gamma^*\).

Illustration of convergence:

The Update Step — MLE Interpretation

Remark (Remark 3.5 in de Boer et al.): the update formula is exactly the maximum likelihood estimator of \(\mathbf{v}\) computed from only the elite samples.

For the Normal family (one component \(j\)):

\[\hat{\mu}_{t,j} = \frac{\sum_{i : G(\mathbf{X}_i) \leq \hat{\gamma}_t} X_{ij}}{\left| E_t \right|}, \qquad \hat{\sigma}_{t,j}^2 = \frac{\sum_{i : G(\mathbf{X}_i) \leq \hat{\gamma}_t} (X_{ij} - \hat{\mu}_{t,j})^2}{\left| E_t \right|}\]

Intuition:

• The means \(\hat{\mu}_{t,j}\) are pulled toward the region of low \(G\) values.
• The standard deviations \(\hat{\sigma}_{t,j}\) shrink as the elites cluster more tightly.
• This is a standard sample mean and sample standard deviation — computed only over the elite subset.

Stopping Criteria

Two independent stopping criteria are used in practice:

Criterion 1 — Distribution convergence (for continuous/Normal sampler):

The distribution has converged when all component standard deviations fall below a threshold based on the coefficient of variation (CV):

\[\hat{\sigma}_{t,j} \leq |\hat{\mu}_{t,j}| \cdot \delta_{\text{CV}} \quad \forall j\]

where \(\delta_{\text{CV}}\) is a small positive value (e.g., 0.01–0.05). This indicates the distribution has degenerated to a near-point mass.

Criterion 2 — No improvement in elite solution:

If the best solution found does not improve over \(d\) consecutive iterations (stagnation stopping), the algorithm halts.

In KSL: both criteria are checked; the algorithm stops when either is satisfied.

Initialization of Parameters

Means \(\mu_j\): initialized to the midpoint of the feasible range for each input dimension \(j\).

Standard deviations \(\sigma_j\): initialized from the input range divided by 4, scaled by the initialVariabilityFactor:

\[\sigma_{0,j} = \frac{\text{range}_j}{4} \times \text{initialVariabilityFactor}\]

Rationale:

• Range/4 approximates \(\sigma\) for a uniform distribution over \([\text{lb}_j, \text{ub}_j]\) — giving broad initial coverage.
• initialVariabilityFactor \(> 1\) increases spread (more exploration); \(< 1\) decreases it (exploitation near the midpoint).

Guard condition: the algorithm rejects an initialization where \(\sigma_{0,j}\) is already below the CV-based convergence threshold. This prevents inadvertently starting in a “converged” state.

Parameter Summary Table

Population Size \(N\) — ceSampleSize

Role: controls the breadth of exploration at each iteration and the statistical quality of the parameter update.

Theoretical guidance (de Boer et al.):

• For problems with \(n\) parameters: \(N \approx c \cdot n\), \(c \in [5, 10]\).
• \(N\) must be large enough that \(N_\text{elite} = \lfloor \rho N \rfloor \geq 2\) (needed to estimate \(\sigma\)).

KSL automatic sizing (recommendCESampleSize):

Uses the Chen-Kelton (1999) formula to find the \(N\) sufficient to estimate the \(\rho\)-quantile at a specified confidence level:

\[N \geq \frac{z_{\alpha/2}^2 \, \rho (1-\rho)}{\epsilon^2}\]

where \(\epsilon\) is the desired half-width on the proportion estimate. Defaults: \(\rho=0.1\), 95% confidence, \(\epsilon=0.1\).

Effect of \(N\): - Too small \(\Rightarrow\) noisy elite set, unstable parameter updates, premature convergence. - Too large \(\Rightarrow\) excessive oracle calls per iteration, slow progress.

Smoother Parameters \(\alpha_\mu\), \(\alpha_\sigma\)

Role: dampen oscillations in the parameter update and prevent premature distribution collapse.

Effect of \(\alpha\) (applies to both mean and SD smoothers):

Practical tip: if the solver converges too quickly to a poor solution, increase \(\alpha\) (more smoothing). If convergence is very slow, decrease \(\alpha\) slightly.

Class Hierarchy

StochasticSolver  (abstract)
│   ├── rnStream, streamProvider       [random stream management]
│   ├── startingPoint()                [feasible starting point]
│   └── mainIteration()  ← abstract
│
└── CrossEntropySolver
        │   ├── ceSampler: CESampler       [the probability model]
        │   ├── ceSampleSize, elitePct     [population parameters]
        │   ├── solutionChecker            [stagnation stopping]
        │   ├── initializeIterations()
        │   ├── mainIteration()            [the CE loop body]
        │   └── isStoppingCriteriaSatisfied()
        │
        └── CESampler  (abstract class)
                │
                └── CENormalSampler
                        ├── myMeans[], myStdDevs[]
                        ├── meanSmoother, sdSmoother, cvThreshold
                        ├── sample()
                        ├── updateParameters()
                        └── hasConverged()

Design Philosophy

Separation of concerns:

Key design benefit: CrossEntropySolver is distribution-agnostic. Any subclass of CESampler (e.g., a Bernoulli sampler for discrete problems) can be plugged in without modifying the solver.

When a sampler is handed to the solver it is attached: adopted onto the solver’s own stream provider on a stream distinct from the solver’s base stream, keeping random number management consistent throughout the hierarchy.

Interface-to-Algorithm Mapping

Note: the base class makes no assumption about which parametric family is used. It only specifies what operations any CE sampler must support.

Constructor and Key Properties

class CENormalSampler(
    problemDefinition: ProblemDefinition,
    meanSmoother: Double  = defaultMeanSmoother,     // α_μ = 0.85
    sdSmoother: Double    = defaultStdDevSmoother,   // α_σ = 0.85
    coefficientOfVariationThreshold: Double          // δ_CV = 0.03
                      = defaultCoefficientOfVariationThreshold,
    streamNum: Int = 0,
    streamProvider: RNStreamProviderIfc = RNStreamProvider(),
) : CESampler(problemDefinition, streamNum, streamProvider)

Internal state:

private val myMeans:   DoubleArray = DoubleArray(dimension) { 1.0 }  // μ_j
private val myStdDevs: DoubleArray = DoubleArray(dimension) { 1.0 }  // σ_j
private val myEliteStats = List(dimension) { Statistic() }           // running stats

All default values are overridable at the companion object level for global effect, or at the instance level for targeted control.

initializeParameters() — Setting \(\hat{\mathbf{v}}_0\)

override fun initializeParameters(values: DoubleArray) {
    // Assign initial means from the supplied starting point (e.g. midpoints)
    for (i in values.indices) {
        require(values[i].isFinite())
        myMeans[i] = values[i]
        myEliteStats[i].reset()
    }
    // Assign initial standard deviations from input ranges
    val ranges = problemDefinition.inputRanges
    for (i in values.indices) {
        myStdDevs[i] = (ranges[i] / 4.0) * initialVariabilityFactor
        val stdDevThreshold = abs(myMeans[i]) * cvThreshold
        require(myStdDevs[i] > stdDevThreshold) {
            "Initial σ (${myStdDevs[i]}) for ${problemDefinition.inputNames[i]}" +
            " already below stopping threshold ($stdDevThreshold)."
        }
    }
}

Algorithm mapping:

• values \(\leftarrow\) starting point supplied by CrossEntropySolver.initializeIterations().
• \(\sigma_{0,j} = \text{range}_j / 4 \times \text{initialVariabilityFactor}\) implements the range-based initialization.
• The require guard prevents a degenerate start.

sample() — Drawing from \(f(\cdot;\hat{\mathbf{v}}_{t-1})\)

override fun sample(array: DoubleArray) {
    require(array.size == dimension)
    for (i in array.indices) {
        // Draw X_j ~ Normal(μ_j, σ_j²)
        // Note: rNormal takes (mean, variance) — σ² not σ
        array[i] = rnStream.rNormal(myMeans[i], myStdDevs[i] * myStdDevs[i])
    }
}

Algorithm mapping: implements Step 1 of Algorithm 3.3 — one component of a single candidate \(\mathbf{X}_i\) is drawn per call. CrossEntropySolver.mainIteration() calls this \(N\) times via the inherited MVSampleIfc.sample(n).

Important: rNormal(mean, variance) takes the variance \(\sigma^2\) as its second argument, not the standard deviation. The squaring myStdDevs[i] * myStdDevs[i] is intentional.

updateParameters() — Computing \(\hat{\mathbf{v}}_t\)

override fun updateParameters(elites: List<DoubleArray>) {
    require(elites.isNotEmpty())
    // Reset and collect statistics over the elite set
    myEliteStats.forEach { it.reset() }
    for (array in elites) {
        for (i in array.indices) {
            myEliteStats[i].collect(array[i])   // accumulate X_ij values
        }
    }
    // Apply smoothed update for means and standard deviations
    for (i in myMeans.indices) {
        // μ_t = α_μ · μ_{t-1} + (1-α_μ) · mean_elite
        myMeans[i] = meanSmoother * myMeans[i] +
                     (1.0 - meanSmoother) * myEliteStats[i].average
        if (myEliteStats[i].count >= 2) {
            // σ_t = α_σ · σ_{t-1} + (1-α_σ) · sd_elite
            myStdDevs[i] = sdSmoother * myStdDevs[i] +
                           (1.0 - sdSmoother) * myEliteStats[i].standardDeviation
        }
        // If count < 2: skip σ update (cannot estimate variance from one point)
    }
}

Algorithm mapping: implements Steps 3–5 of Algorithm 3.3. The elites list is the output of CrossEntropySolver.findEliteSolutions().

hasConverged() — CV-Based Stopping

override fun hasConverged(): Boolean {
    for (i in myStdDevs.indices) {
        // If any dimension's σ exceeds its threshold, NOT converged
        if (myStdDevs[i] > abs(myMeans[i]) * cvThreshold) return false
    }
    return true  // All dimensions satisfied: distribution has collapsed
}

// Inspection property: shows threshold for each dimension
val stdDeviationThresholds: DoubleArray
    get() = DoubleArray(dimension) { abs(myMeans[it]) * cvThreshold }

Algorithm mapping: implements the distribution-convergence branch of Step 6 of Algorithm 3.3.

CV condition per dimension \(j\):

\[\hat{\sigma}_{t,j} \leq |\hat{\mu}_{t,j}| \cdot \delta_\text{CV}\]

This is equivalent to requiring the coefficient of variation \(\text{CV}_j = \hat{\sigma}_{t,j} / |\hat{\mu}_{t,j}|\) to be below \(\delta_\text{CV}\) for every dimension simultaneously.

toString() as a Diagnostic

override fun toString(): String {
    return buildString {
        appendLine("CENormalSampler")
        appendLine("streamNumber = $streamNumber")
        appendLine("dimension = $dimension")
        appendLine("variabilityFactor = $initialVariabilityFactor")
        appendLine("meanSmoother = $meanSmoother")
        appendLine("sdSmoother = $sdSmoother")
        appendLine("coefficient of variation threshold = $cvThreshold")
        appendLine("mean values = ${myMeans.contentToString()}")
        appendLine("standard deviations = ${myStdDevs.contentToString()}")
        appendLine("std deviation thresholds = ${stdDeviationThresholds.contentToString()}")
        append("Have standard deviation thresholds converge? = ${hasConverged()}")
    }
}

Printing the sampler mid-run shows current \(\boldsymbol{\mu}\), \(\boldsymbol{\sigma}\), and their convergence thresholds — useful for diagnosing premature or stalled convergence.

initializeIterations() — Algorithm Initialization

override fun initializeIterations() {
    // Obtain the starting point (user-supplied or randomly generated)
    val initialPoint = startingPoint ?: startingPoint()

    // Initialize the CE sampler's parameters from the starting point
    // This sets μ_0 = starting input values, σ_0 from input ranges
    ceSampler.initializeParameters(initialPoint.inputValues)

    myEliteSolutions.clear()
    solutionChecker.clear()
}

Algorithm mapping: implements the initialization step of Algorithm 3.3 — sets \(\hat{\mathbf{v}}_0 = \mathbf{u}\).

The starting point is either: - User-supplied via startingPoint property, or - Generated by startingPoint() which uses StochasticSolver’s feasibility-sampling mechanism.

mainIteration() — The CE Loop Body (Part 1)

override fun mainIteration() {
    // ── Step 1: Generate the CE population ──────────────────────────────
    // Draws N candidates X_1,...,X_N ~ f(·; v_{t-1})
    val points = ceSampler.sample(sampleSize())         // N × d matrix

    // Convert raw Double arrays to InputMap objects (respects feasibility)
    val inputs = problemDefinition.convertPointsToInputs(points)

    // ── Step 2: Evaluate each candidate via the simulation oracle ───────
    // Each evaluation uses `replicationsPerEvaluation` replications
    val evaluations = requestEvaluations(inputs)
    val results = evaluations.values.toList()

    if (results.isEmpty()) return  // Safety: no evaluations returned

Algorithm mapping: Steps 1–2 of Algorithm 3.3. The oracle evaluations are the dominant computational cost.

mainIteration() — The CE Loop Body (Part 2)

    // ── Step 3: Identify elite set ───────────────────────────────────────
    // Sort by performance; retain top eliteSize() solutions
    // eliteSize() = max(⌈elitePct · ceSampleSize⌉, defaultMinEliteSize)
    myEliteSolutions = findEliteSolutions(results)      // sorted, top-ρN

    // ── Steps 4–5: Update sampler parameters (with smoothing) ───────────
    // Extract the input coordinate arrays from elite Solution objects
    val elitePoints = extractSolutionInputPoints(myEliteSolutions)

    // Delegates to CENormalSampler.updateParameters():
    //   - computes mean and SD of elitePoints per dimension
    //   - applies exponential smoothing
    ceSampler.updateParameters(elitePoints)

    // ── Track best solution ──────────────────────────────────────────────
    currentSolution = myEliteSolutions.first()   // best among elites
    solutionChecker.captureSolution(currentSolution)
}

Algorithm mapping: Steps 3–5 and best-solution tracking of Algorithm 3.3.

isStoppingCriteriaSatisfied() — Step 6

override fun isStoppingCriteriaSatisfied(): Boolean {
    return solutionQualityEvaluator
               ?.isStoppingCriteriaReached(this)
           ?: checkForConvergence()
}

private fun checkForConvergence(): Boolean {
    // Criterion 1: stagnation — no improvement in best solution
    //              for `noImproveThreshold` consecutive iterations
    // Criterion 2: distribution convergence — CENormalSampler.hasConverged()
    //              all σ_j ≤ |μ_j| · δ_CV
    return solutionChecker.checkSolutions() || ceSampler.hasConverged()
}

Two independent stopping criteria — either triggers termination:

Elite Size Determination — eliteSize()

fun eliteSize(): Int {
    // If a custom function is provided, use it (allows dynamic elite sizing)
    return eliteSizeFn?.eliteSize(this)
        // Otherwise: max(⌈ρ·N⌉, defaultMinEliteSize)
        ?: maxOf(ceil(elitePct * ceSampleSize).toInt(), defaultMinEliteSize)
}

The EliteSizeIfc and SampleSizeFnIfc functional interfaces allow users to supply dynamic sizing strategies:

// Example: double the elite size after iteration 10
solver.eliteSizeFn = EliteSizeIfc { ceSolver ->
    if (ceSolver.iterationCounter > 10)
        2 * ceil(ceSolver.elitePct * ceSolver.ceSampleSize).toInt()
    else
        ceil(ceSolver.elitePct * ceSolver.ceSampleSize).toInt()
}

This supports adaptive variants analogous to the FACE algorithm described in de Boer et al.

Basic Constructor Call

// Create a problem and evaluator (problem-specific)
val problem   = ProblemDefinition(...)
val evaluator = MySimulationEvaluator(...)

// Build a CENormalSampler — tune algorithmic parameters here
val sampler = CENormalSampler(
    problemDefinition              = problem,
    meanSmoother                   = 0.85,   // α_μ: smoothing on means
    sdSmoother                     = 0.85,   // α_σ: smoothing on SDs
    coefficientOfVariationThreshold = 0.03   // δ_CV: convergence threshold
)
sampler.initialVariabilityFactor = 1.0       // starting σ multiplier

// Build the solver
val solver = CrossEntropySolver(
    problemDefinition        = problem,
    evaluator                = evaluator,
    ceSampler                = sampler,
    maxIterations            = 100,          // hard iteration limit
    replicationsPerEvaluation = 10           // oracle noise level
)

// Tune population parameters
solver.ceSampleSize = 50                     // N
solver.elitePct     = 0.10                   // ρ = 10%

solver.runAllIterations()
println(solver.bestSolution)

Global Defaults via Companion Objects

Both CENormalSampler and CrossEntropySolver expose companion object properties that act as global defaults for all subsequently created instances:

// Change defaults before creating any solvers
CENormalSampler.defaultMeanSmoother                   = 0.80
CENormalSampler.defaultStdDevSmoother                 = 0.80
CENormalSampler.defaultCoefficientOfVariationThreshold = 0.02
CENormalSampler.defaultInitialVariabilityFactor        = 1.2

CrossEntropySolver.defaultElitePct                    = 0.05
CrossEntropySolver.defaultMinEliteSize                = 10
CrossEntropySolver.defaultNoImproveThresholdForCE     = 7
CrossEntropySolver.defaultMaxCESampleSize             = 200
CrossEntropySolver.ceDefaultMaxIterations             = 150

Use case: when running a batch of experiments with the same parameter profile, set the companion object defaults once rather than configuring each solver instance individually.

Diagnostic Guidance

toString() for Full Solver Diagnostics

println(solver.toString())
// CrossEntropySolver(
//     elitePct = 0.1,
//     ceSampleSize = 50,
//     eliteSize = 5,
//     eliteSizeFn = None,
//     sampleSizeFn = None,
//     noImproveThreshold = 5,
//     ceSampler = CENormalSampler,
//     base = StochasticSolver(
//         streamNumber = 1, ...
//     )
// )

println(solver.ceSampler.toString())
// CENormalSampler
// streamNumber = 1
// dimension = 3
// variabilityFactor = 1.0
// meanSmoother = 0.85
// sdSmoother = 0.85
// coefficient of variation threshold = 0.03
// mean values = [4.12, 7.83, 2.95]
// standard deviations = [0.21, 0.18, 0.09]
// std deviation thresholds = [0.12, 0.24, 0.09]
// Have standard deviation thresholds converge? = false

Use these outputs at the end of each run — or after early termination — to understand what the algorithm found and how close to convergence it was.

CE Algorithm \(\leftrightarrow\) KSL Implementation

Key Takeaways

1. CE reformulates optimization as adaptive importance sampling. The sampling distribution \(f(\cdot;\hat{\mathbf{v}}_t)\) is steered toward the optimum by repeatedly fitting it to the elite performers.

2. The update step is analytic for the Normal family. Means and standard deviations of the elite set are computed as sample statistics — no numerical optimization required inside the loop.

3. Smoothed updating is essential. Without it, the distribution collapses prematurely. The default smoothers (\(\alpha = 0.85\)) represent a good starting point for most problems.

4. Two stopping criteria work together. CV-based convergence confirms the distribution has degenerated to a solution; stagnation-based stopping provides a safety net when convergence is slow.

5. CrossEntropySolver is distribution-agnostic. Swapping CENormalSampler for any other CESampler subclass immediately generalizes CE to discrete, integer, or custom domains.

6. Diagnosis starts with toString(). Printing the sampler state reveals current means, standard deviations, and thresholds — enabling rapid identification of parameter tuning needs.

When to Choose CE

References

• de Boer, P.-T., Kroese, D. P., Mannor, S., and Rubinstein, R. Y. (2005). A Tutorial on the Cross-Entropy Method. Annals of Operations Research, 134(1), 19–67.

• Rubinstein, R. Y. (1999). The simulated entropy method for combinatorial and continuous optimization. Methodology and Computing in Applied Probability, 2, 127–190.

• Rubinstein, R. Y. and Kroese, D. P. (2004). The Cross-Entropy Method: A Unified Approach to Combinatorial Optimization, Monte-Carlo Simulation and Machine Learning. Springer.

• Chen, H. and Kelton, W. D. (1999). Quantile and tolerance-interval estimation in simulation. Proceedings of the 1999 Winter Simulation Conference.

• Rossetti, M. D. KSL — Kotlin Simulation Library. ksl.simopt.solvers.algorithms.