Integer-Ordered Simulation Optimization

The General SO Problem:

\[\min_{x \in X \subseteq \mathbb{Z}^d} \; g(x)\]

where:

• \(g : x \mapsto \mathbb{R}\) is a black-box objective, observable only via a stochastic simulation oracle.
• At any \(x\), only a noisy estimate \(\hat{g}_m(x)\) is available from \(m\) replications.
• As \(m \to \infty\): \(\hat{g}_m(x) \to g(x)\) with probability 1.
• The decision variables \(x = (x_1, x_2, \ldots, x_d)\) are integer-valued.

Practical examples: inventory reorder points and quantities, number of machines/servers, staffing levels, buffer sizes.

Why Is This Hard?

Three compounding challenges:

1. Stochastic noise — we never observe \(g(x)\) exactly; only a noisy estimate \(\hat{g}_m(x)\).
2. No analytical gradient — \(g\) is implicit in a simulation; derivatives are not available in closed form.
3. Integer domain — continuous gradient methods cannot be applied directly; combinatorial search is expensive.

The key insight of R-SPLINE: construct a continuous relaxation of \(\hat{g}_m(\cdot)\) over the integer lattice using piecewise-linear interpolation, apply gradient search on that relaxation, then snap back to integers.

Classification of SO Algorithms

• R-SPLINE targets local optima on integer-ordered, potentially infinite feasible spaces.
• Of 50 benchmark SO problems in the literature, 25 fall in exactly this category (Wang et al. 2013).
• Primary competitors: COMPASS, ISC (Industrial Strength COMPASS), ISC-AHA.

What Does R-SPLINE Stand For?

R-SPLINE = Retrospective Search with Piecewise-Linear Interpolation and Neighborhood Enumeration

Defining Neighborhoods

Definition (Neighborhood): The neighborhood \(N(\mathbf{0})\) of the \(d\)-dimensional origin is any subset of \(\mathbb{Z}^d\) containing the origin. The neighborhood of any \(x \in \mathbb{Z}^d\) is: \[N(x) = \{y : (y - x) \in N(\mathbf{0})\}\]

Two important special cases used in R-SPLINE:

• \(N_1(x)\): the \(2d\) integer points exactly one unit away from \(x\) (von Neumann neighborhood). \[N_1(x) = \{x \pm e_i : i = 1, \ldots, d\}\]
• \(N_2(x)\): all integer points within 2 units of \(x\).

R-SPLINE targets \(N_1\)-local minima — the smallest computationally tractable neighborhood.

N1-Local Minimum

Definition: A point \(x^* \in = X\) is an \(N_1\)-local minimizer of \(g\) if: \[g(x) \geq g(x^*) \quad \forall x \in N_1(x^*) \cap X\]

That is, \(x^*\) is no worse than any of its \(2d\) immediate neighbors.

Key properties:

• The set of \(N_1\)-local minima \(M^*(N_1)\) contains all global minima.
• Larger neighborhoods \(\Rightarrow\) fewer local minima (higher quality but costlier to verify).
• R-SPLINE uses \(N_1\) because NE requires evaluating all \(2d\) neighbors — manageable even in high \(d\).
• SPLI compensates for the small neighborhood by allowing large, off-axis moves.

In KSL: The Von Neumann Neighborhood

The KSL implementation uses a vonNeumannNeighborhoodFinder by default, which directly realizes the \(N_1\) neighborhood:

var neighborhoodFinder: NeighborhoodFinderIfc =
    problemDefinition.vonNeumannNeighborhoodFinder()

• The neighborhood finder returns the set of all input maps one unit away in each coordinate direction.
• Only input-feasible neighbors (satisfying bounds and deterministic constraints) are retained.
• Users can substitute a different neighborhood finder if a larger structure is preferred.

R-SPLINE Architecture Overview

R-SPLINE (outer loop — Algorithm 1)
│
│  for k = 1, 2, 3, ...
│    solve sample-path problem P_k using SPLINE
│    warm-start next iteration with X*_k
│
└─► SPLINE (Algorithm 2)
     │
     │  repeat until N1-local optimum found or budget b_k exhausted:
     │    SPLI → NE → SPLI → NE → ...
     │
     ├─► SPLI (Algorithm 4)   — gradient line search on continuous relaxation
     │     └─► PLI (Algorithm 3) — piecewise-linear interpolation + gradient
     │           └─► PERTURB      — shifts integer point into a simplex interior
     │
     └─► NE                   — enumerate all N1 neighbors, pick the best

Each level has a distinct role. Together they combine the efficiency of gradient search with the correctness guarantee of discrete local-optimality checking.

The Retrospective Framework (Algorithm 1)

Core idea: solve a sequence of sample-path problems \(P_k\), each with a larger sample size \(m_k\):

\[(P_k): \text{ Find } X^*_k \in X \text{ such that } \hat{g}_{m_k}(x) \geq \hat{g}_{m_k}(X^*_k) \quad \forall x \in N_1(X^*_k)\]

Algorithm 1 pseudocode:

Set X*_0 = x0  (initial solution)
for k = 1, 2, 3, ...
    X*_k = SPLINE(X*_{k-1}, m_k, b_k)
end for

Three distinguishing features:

1. Independent sample paths — different random numbers used each iteration \(k\).
2. Increasing sample sizes — \(m_k \uparrow \infty\), so each \(P_k\) is a better approximation of \(g\).
3. Warm starts — \(X^*_{k-1}\) seeds iteration \(k\), making later iterations fast.

Retrospective Framework: Why It Works

Early iterations (\(k\) small, \(m_k\) small):

• Computation per point is cheap.
• SPLINE explores broadly; many oracle calls acceptable because each is inexpensive.
• Purpose: quickly find a good neighborhood of the true local minimum.

Later iterations (\(k\) large, \(m_k\) large):

• Warm start \(X^*_{k-1}\) is already near the true local minimum.
• SPLINE needs few oracle calls (sample path agrees with \(g\) closely).
• Purpose: confirm the local minimum with high statistical confidence.

Efficiency gain: most of the computation budget is spent in early iterations where evaluations are cheap.

Mapping to KSL: mainIteration()

Algorithm 1’s loop body maps directly to mainIteration() in RSplineSolver:

override fun mainIteration() {
    // Call SPLINE for the next solution using
    // current sample size (m_k) and SPLINE call limit (b_k).
    val splineSolution = spline(
        currentSolution,
        rSPLINESampleSize,   // this is m_k
        splineCallLimit      // this is b_k
    )
    currentSolution = splineSolution
    solutionChecker.captureSolution(currentSolution)
}

• rSPLINESampleSize returns \(m_k\) from the FixedGrowthRateReplicationSchedule.
• splineCallLimit computes \(b_k = b_0 \cdot (1 + r_b)^{k-1}\), capped at maxSplineCallLimit.
• solutionChecker accumulates solutions for the KSL-specific convergence stopping criterion.

The Sample Size Schedule

The paper specifies a 10% static growth rate as the default: \[m_{k+1} = \lceil 1.1 \cdot m_k \rceil\]

In KSL, this is implemented by FixedGrowthRateReplicationSchedule:

constructor(
    initialNumReps: Int = defaultInitialSampleSize,      // m_1, default = 8
    sampleSizeGrowthRate: Double = defaultReplicationGrowthRate,  // 0.10
    maxNumReplications: Int = defaultMaxNumReplications  // upper cap
)

The schedule satisfies the theoretical requirement that \(m_k \uparrow \infty\) faster than logarithmically, ensuring convergence. The growth rate is very mild — only 10% more replications per iteration.

SPLINE: The Sample-Path Solver (Algorithm 2)

For a fixed \(k\), SPLINE finds an \(N_1\)-local minimum of the sample-path objective \(\hat{g}_{m_k}(\cdot)\):

Inputs: initial solution X_new, sample size m_k, oracle call budget b_k
Set oracle call counter n = 0
repeat
    [n', X_old] = SPLI(X_new, m_k)     // gradient line search
    if g_hat(X_old) > g_hat(X_new):
        X_old = X_new                   // SPLI cannot cause harm
    [n'', X_new] = NE(X_old, m_k)      // neighborhood check
    n = n + n' + n''
until X_new == X_old  OR  n > b_k

Termination: SPLINE stops when NE confirms \(X_{old}\) is \(N_1\)-locally optimal (no neighbor is better), or the call budget \(b_k\) is exhausted.

SPLINE in KSL: spline() Function

private fun spline(initSolution: Solution,
                   sampleSize: Int, splineCallLimit: Int): Solution {
    var newSolution = /* re-evaluate at m_k if needed */
    val startingSolution = newSolution
    for (i in 1..splineCallLimit) {
        val spliResults = searchPiecewiseLinearInterpolation(newSolution, sampleSize)
        // SPLI cannot cause harm: keep current if SPLI didn't improve
        val neStartingSolution = if (compare(spliResults.solution, newSolution) <= 0)
            spliResults.solution else newSolution
        val neSearchResults = neighborhoodSearch(neStartingSolution, sampleSize)
        newSolution = neSearchResults.solution
        // If NE confirmed local optimality: break
        if (neStartingSolution == neSearchResults.solution) break
    }
    // Return best feasible solution found; fall back to startingSolution if no improvement
    return if (compare(startingSolution, newSolution) < 0) startingSolution
           else if (newSolution.isInputFeasible()) newSolution
           else startingSolution
}

Key implementation notes:

• The loop bound splineCallLimit corresponds to \(b_k\) in the paper. Safety check: if SPLINE returns an infeasible solution, the starting solution is returned. The re-evaluation step ensures the initial solution uses the current \(m_k\) sample size.

Neighborhood Enumeration: NE

NE is the simplest component — it exhaustively evaluates all \(N_1\) neighbors:

Evaluate all feasible x in N1(X_old) using m_k replications each.
Case 1: Unique best neighbor x better than X_old → return x.
Case 2: X_old tied for best → return X_old (local optimum confirmed).
Case 3: Multiple neighbors tied for best and better → return any one.

In KSL (neighborhoodSearch()):

val neighborhood = neighborhoodFinder.neighborhood(solution.inputMap, this)
val feasible = neighborhood.filter { it.isInputFeasible() }.toSet()
val evaluations = requestEvaluations(feasible, sampleSize)
val candidate = results.minOf { it }
return if (compare(solution, candidate) < 0)
    NESearchResults(results.size * sampleSize, solution)  // already local min
else
    NESearchResults(results.size * sampleSize, candidate) // improved

NE is the local optimality verifier — SPLI does the exploration, NE does the confirmation.

Why Interpolate?

The objective \(\hat{g}_{m_k}(\cdot)\) is only defined on integer points \(X \subseteq \mathbb{Z}^d\).

To perform a gradient-based line search, we need a continuous extension \(g_{m_k}(\cdot) : \mathbb{R}^d \to \mathbb{R}\) that:

1. Agrees with \(\hat{g}_{m_k}(\cdot)\) at every integer point.
2. Is continuous everywhere (so the gradient is well-defined almost everywhere).
3. Is cheap to compute (only \(d+1\) oracle calls per simplex evaluation).
4. Has local solutions that encompass the \(N_1\)-optimal integer solutions.

The chosen method: simplicial piecewise-linear interpolation. Every non-integer point \(x \in \mathbb{R}^d \setminus \mathbb{Z}^d\) lies inside a simplex formed by \(d+1\) nearby integer vertices.

The Simplex Construction

Step 1: Find the floor vertex \(x^0\): \[x^0_i = \lfloor x_i \rfloor, \quad i = 1, \ldots, d\]

Step 2: Compute fractional parts \(z_i = x_i - x^0_i \in [0,1)\).

Step 3: Sort the fractional parts in descending order. Let \(p(1), p(2), \ldots, p(d)\) be the permutation such that: \[1 = z_{p(0)} > z_{p(1)} \geq z_{p(2)} \geq \cdots \geq z_{p(d)} \geq z_{p(d+1)} = 0\]

Step 4: Build the remaining \(d\) vertices by sequentially adding unit vectors: \[x^i = x^{i-1} + e_{p(i)}, \quad i = 1, 2, \ldots, d\]

Step 5: Compute vertex weights (convex combination coefficients): \[w_i = z_{p(i)} - z_{p(i+1)}, \quad i = 0, 1, \ldots, d\]

Note: \(\sum_{i=0}^d w_i = z_{p(0)} - z_{p(d+1)} = 1 - 0 = 1\). ✓

Worked Example: \(x = (1.8, 2.3, 3.6)\)

Step 1: Floor vertex: \(x^0 = (1, 2, 3)\).

Step 2: Fractional parts: \(z = (0.8, 0.3, 0.6)\).

Step 3: Sort descending: \(z_{p(1)} = 0.8\) (dim 1), \(z_{p(2)} = 0.6\) (dim 3), \(z_{p(3)} = 0.3\) (dim 2). Sorted sequence: \(1.0, \; 0.8, \; 0.6, \; 0.3, \; 0.0\)

Step 4: Vertices: \[x^0 = (1,2,3), \quad x^1 = (2,2,3), \quad x^2 = (2,2,4), \quad x^3 = (2,3,4)\]

Step 5: Weights: \(w_0 = 1.0 - 0.8 = 0.2\), \(w_1 = 0.8 - 0.6 = 0.2\), \(w_2 = 0.6 - 0.3 = 0.3\), \(w_3 = 0.3 - 0.0 = 0.3\)

Interpolated value: \[g_{m_k}(1.8, 2.3, 3.6) = 0.2 \cdot \hat{g}(1,2,3) + 0.2 \cdot \hat{g}(2,2,3) + 0.3 \cdot \hat{g}(2,2,4) + 0.3 \cdot \hat{g}(2,3,4)\]

Gradient Estimation from the Simplex

Once all \(d+1\) vertices are evaluated, the PLI provides a phantom gradient estimate:

\[\hat{\gamma}_{p(i)} = \hat{g}_{m_k}(x^i) - \hat{g}_{m_k}(x^{i-1}), \quad i = 1, 2, \ldots, d\]

Interpretation: Each component \(\hat{\gamma}_{p(i)}\) measures the change in the estimated objective when moving one unit in the \(p(i)\)-th coordinate direction, at the current simplex location.

Properties of the gradient estimate:

• Available only when all \(d+1\) vertices are feasible.
• If any vertex is infeasible, gradient computation is skipped (SPLI falls back gracefully).
• Uses common random numbers across all simplex vertices → reduces variance of comparisons.

PLI: Algorithm 3 Summary

\[g_{m_k}(x) = \frac{\sum_{i=0}^{d} w_i \cdot \hat{g}_{m_k}(x^i) \cdot \mathbf{1}[x^i \text{ feasible}]}{\sum_{i=0}^{d} w_i \cdot \mathbf{1}[x^i \text{ feasible}]}\]

Two key properties of \(g_{m_k}(\cdot)\):

1. Consistency: \(g_{m_k}(x) = \hat{g}_{m_k}(x)\) for every \(x \in X\) (agrees with oracle at integer points).
2. Encompassing: Every \(N_1\)-local minimum of \(\hat{g}_{m_k}(\cdot)\) is also a local minimum of \(g_{m_k}(\cdot)\).

Property 2 is why we can search for a continuous local minimum and know we have not missed any discrete local minima.

piecewiseLinearSimplex() in KSL

This companion object function implements the simplex construction (Algorithm 3, Steps 1–5):

fun piecewiseLinearSimplex(point: DoubleArray): SimplexData {
    // Step 1: floor vertex x0
    val x0 = DoubleArray(point.size) { floor(point[it]) }
    vertices.add(x0)
    // Step 2: fractional parts z
    val z = DoubleArray(point.size) { point[it] - x0[it] }
    // Step 3: sort indices by descending fractional part
    val zSortedIndices = z.sortIndices(descending = true)
    // Step 4: build remaining vertices
    var np = x0
    for ((index, _) in zSortedIndices.withIndex()) {
        val ei = DoubleArray(z.size); ei[zSortedIndices[index]] = 1.0
        val nx = KSLArrays.addElements(np, ei)
        np = nx; vertices.add(nx)
    }
    // Step 5: compute weights w_i = z_p(i) - z_p(i+1)
    val zList = mutableListOf(1.0); zSortedIndices.forEach { zList.add(z[it]) }; zList.add(0.0)
    val w = DoubleArray(z.size + 1) { zList[it] - zList[it + 1] }
    return SimplexData(point, z, zSortedIndices, zList, vertices, w)
}

Note: This is a pure function — no simulation calls occur here. It is static (companion object), making it easily testable in isolation.

piecewiseLinearInterpolation() in KSL

The piecewiseLinearInterpolation() private function orchestrates the full PLI (Algorithm 3 Steps 6–13):

private fun piecewiseLinearInterpolation(solution, sampleSize, perturbation): PLIResults {
    val x = solution.inputMap.inputValues     // current integer point
    val point = perturb(x, perturbation)      // shift into simplex interior (PERTURB)
    val simplexData = piecewiseLinearSimplex(point) // geometry (Algorithm 3 steps 1-5)
    val feasibleInputs = filterToFeasibleInputs(simplexData.simplexPoints)
    // Evaluate feasible vertices using CRN
    val evaluations = requestEvaluationsWithCRN(feasibleInputs.keys, sampleSize)
    // Gradient computation (Algorithm 3 step 12)
    val gradients = DoubleArray(d)
    for ((i, indexValue) in simplexData.sortedFractionIndices.withIndex()) {
        gradients[indexValue] = results[i+1].penalizedObjFncValue
                              - results[i].penalizedObjFncValue
    }
    return PLIResults(numOracleCalls, gradients, bestSolution)
}

CRN note: The requestEvaluationsWithCRN() call synchronizes random number streams across all simplex vertices — directly implementing the paper’s common random numbers requirement.

The PERTURB Step

Why perturb? The PLI algorithm requires a non-integer input point \(x_1 \in \mathbb{R}^d \setminus \mathbb{Z}^d\) to define a simplex. The current best solution \(X_{best}\) is always an integer point, so it must be randomly displaced.

The PERTURB function adds a small uniform random perturbation to each coordinate:

fun addRandomPerturbation(point: DoubleArray,
                          perturbationFactor: Double,
                          rnStream: RNStreamIfc): DoubleArray {
    for (i in point.indices) {
        point[i] = point[i] + rnStream.rUniform(-perturbationFactor, perturbationFactor)
    }
    return point
}

• perturbationFactor (default 0.15) must be in \((0, 1)\) to keep the perturbed point strictly within the surrounding hypercube.
• Perturbation randomizes which simplex \(X_{best}\) is assigned to, providing variation in the gradient direction across SPLI iterations.

SPLI Overview (Algorithm 4)

SPLI performs a gradient-based line search on the continuous relaxation \(g_{m_k}(\cdot)\). It returns \(X_{best}\): the best integer solution found.

Set X_best = x_0,  n' = 0
loop (new line search with new gradient estimate):
    x_1 = PERTURB(X_best)
    Call PLI(x_1, m_k) → get g_mk(x_1) and gradient γ̂
    If ||γ̂|| = 0 or undefined: STOP (no direction available)
    i = 0;  x_0 = X_best
    loop (line search steps with doubling step size):
        i++
        s = c^(i-1) * s_0;   x_1 = x_0 - s * (γ̂ / ||γ̂||)
        Round x_1 to nearest integer.
        If x_1 infeasible: STOP
        Simulate at x_1; update n'
        If i ≤ 2 and no improvement: STOP (wrong direction)
        If i > 2 and no improvement: BREAK (continue outer loop)
    go back to top (new gradient estimate)

SPLI: Step Size Doubling

The line search uses an exponentially increasing step size: \[s_i = s_0 \cdot c^{i-1}, \quad s_0 = 2.0, \; c = 2.0\]

• \(s_0 = 2.0\) ensures the first step moves outside the immediate hypercube (avoids returning to \(x^0\) trivially).
• Doubling allows large moves — directly compensating for the small \(N_1\) neighborhood.
• Movement direction: negative gradient (steepest descent): \(x_1 = x_0 - s \cdot \hat{\gamma} / \|\hat{\gamma}\|\)

SPLI: Four Stopping Conditions

SPLI halts a line search early if any of the following occur:

After step \(i > 2\), a non-improving step breaks the inner loop to start a new gradient estimate rather than stopping entirely — allowing SPLI to recover and search in a different direction.

SPLI in KSL: searchPiecewiseLinearInterpolation()

private fun searchPiecewiseLinearInterpolation(solution, sampleSize): SPLIResults {
    var bestSoln = solution
    for (j in 1..spliMaxIterations) {           // outer: new gradient estimates
        val pliResults = piecewiseLinearInterpolation(bestSoln, sampleSize, perturbation)
        bestSoln = /* update if PLI improved */
        if (pliResults.gradients == null) return SPLIResults(numOracleCalls, bestSoln)
        val direction = pliResults.gradients.direction()
        val x0 = bestSoln.inputMap.inputValues
        for (i in 1..lineSearchIterMax) {        // inner: doubling line search
            val stepSize = s0 * c.pow(i - 1)
            val x1 = x0 - stepSize * direction
            val inputs = problemDefinition.toInputMap(x1) // rounds to integer
            if (!inputs.isInputFeasible()) return SPLIResults(numOracleCalls, bestSoln)
            val x1Solution = requestEvaluation(inputs, sampleSize)
            if (i <= 2 && compare(x1Solution, bestSoln) >= 0)
                return SPLIResults(numOracleCalls, bestSoln)
            if (i >  2 && compare(x1Solution, bestSoln) >= 0) break
            bestSoln = x1Solution
        }
    }
    return SPLIResults(numOracleCalls, bestSoln)
}

Paper vs. KSL: Key Implementation Differences

The KSL implementation adds explicit iteration limits on inner loops to prevent runaway computation in edge cases, while preserving the algorithm’s essential structure.

Convergence Setup

Two assumptions underpin R-SPLINE’s convergence guarantees:

Assumption 1 (Bounded level sets): For every \(x \in X\), there exists \(\delta > 0\) such that the level set \[S(x, \delta) = \{x' \in X : g(x') \leq g(x) + \delta\}\] is finite. This prevents the iterates from “escaping to infinity.”

Assumption 2 (Fast enough probability decay): For \(x \in X\) and \(y \notin S(x,\delta)\), there exists a summable sequence \(\{q_k\}\) such that \(P\{y \in S_k(x)\} \leq q_k\) and \(\sum_{k=1}^\infty q_k < \infty\).

This is a Borel-Cantelli condition ensuring that bad sample-path outcomes become rare fast enough.

When Do the Assumptions Hold?

Lemma 5.1 (Large Deviations): If \(\{\hat{g}_{m_k}(x) - g(x)\}\) satisfies a large-deviation principle with a non-degenerate rate function \(I_x(s)\), and the sample sizes grow faster than logarithmically: \[\limsup_{k \to \infty} m_k^{-1} (\log k)^{1+\varepsilon} = 0 \text{ for some } \varepsilon > 0\] then Assumption 2 holds.

Practical implication: The 10% growth schedule \(m_{k+1} = \lceil 1.1 m_k \rceil\) easily satisfies this (geometric growth \(\gg\) logarithmic). Even \(m_k = (\log k)^2\) is sufficient.

Lemma 5.2 (CLT-based): If a CLT holds for \(\hat{g}_{m_k}\) with uniformly bounded variance, Assumption 2 holds for sample sizes growing faster than quadratically in \(k\).

Main Convergence Results

Theorem 5.4 (Almost-sure convergence):

Under Assumptions 1 and 2, for any starting point \(x_0\): \[P\{X^*_k \notin M^*(N_1) \text{ infinitely often}\} = 0\]

That is, after a finite random number of iterations, all subsequent iterates are true \(N_1\)-local minima, with probability one.

Theorem 5.5 (Exponential convergence rate):

If \(X\) is finite and \(\{\hat{g}_{m_k}(x) - g(x)\}\) satisfies a large-deviation principle:

1. \(P\{M^*_k(N_1) \neq M^*(N_1)\} = O(\exp\{-\beta m_k\})\) for some \(\beta > 0\).

2. If \(m_k\) grows at least linearly: \(P\{X^*_k \notin M^*(N_1)\} = O(\exp\{-\beta k\})\).

The probability of returning a false local minimum decays exponentially in the sample size.

Interpretation of Convergence

Important caveat: R-SPLINE converges to a local minimum — not necessarily the global minimum. When multiple local minima exist with differing objective values, ISC/ISC-AHA (with their global Phase I search) may return better-quality solutions, at the cost of more oracle calls.

Sample Size Schedule Parameters

These parameters control the \(m_k\) sequence — the most theoretically important parameters:

Effect: Larger \(m_1\) gives a more reliable starting point but delays the warm-start advantage. Higher growth rate converges faster theoretically but spends more budget per iteration.

SPLINE Call Budget Parameters

These parameters control \(b_k\) — the per-iteration budget for SPLINE iterations:

Computed as: \[b_k = \min\left(b_{max}, \left\lceil b_0 \cdot (1 + r_b)^{k-1} \right\rceil\right)\]

A growing \(b_k\) is needed for convergence when \(X\) is unbounded — SPLINE must be allowed to search longer as iterations proceed.

SPLI Search Parameters

These control the inner gradient search within each SPLINE iteration:

Advice: The defaults are well-tested. Increase spliMaxIterations or lineSearchIterMax for high-dimensional problems where a single gradient estimate may not provide enough search direction.

Stopping and Outer Loop Parameters

KSL-specific stopping criterion: The paper does not specify when to stop the outer loop. KSL adds a practical stopping rule: if the last \(N\) consecutive solutions are statistically indistinguishable (same inputs AND no significant difference in penalized objective function via confidence interval test), the algorithm terminates early.

Parameter Summary Table

Constructor and Setup

RSplineSolver requires an integer-ordered problem definition. Two constructor forms:

// Convenience constructor — specify growth schedule by its parameters directly
val solver = RSplineSolver(
    problemDefinition = myProblem,
    evaluator = myEvaluator,
    maxIterations = 200,
    initialNumReps = 10,                    // m_1
    sampleSizeGrowthRate = 0.10,            // 10% growth
    maxNumReplications = 500,
    solutionEqualityChecker = InputsAndConfidenceIntervalEquality()
)

// Tune SPLI parameters
solver.perturbation = 0.15
solver.initialLineSearchStepSize = 2.0
solver.lineSearchStepSizeMultiplier = 2.0
solver.spliMaxIterations = 5
solver.lineSearchIterMax = 10

// Run
solver.solve()

Requirement: problemDefinition.isIntegerOrdered must be true — enforced at construction time.

Common Random Numbers in KSL

The paper emphasizes that all evaluations within iteration \(k\) use the same random seed \(\xi_k\), making comparisons between points more reliable.

In KSL, this is handled automatically:

• requestEvaluationsWithCRN() is called for simplex vertex evaluation inside piecewiseLinearInterpolation().
• The rnStream property provides a single random number stream for perturbation.
• The evaluator infrastructure manages CRN alignment across oracle calls within an iteration.

Effect: Common random numbers create positive correlation between \(\hat{g}_{m_k}(x_1)\) and \(\hat{g}_{m_k}(x_2)\) for different \(x_1, x_2\), which increases the probability that comparisons correctly identify the better point.

When to Choose R-SPLINE

R-SPLINE at a Glance

Key Takeaways

1. R-SPLINE solves integer-ordered local SO problems by alternating a continuous gradient search (SPLI) with a discrete local-optimality check (NE), inside a retrospective framework.

2. Piecewise-linear interpolation creates a continuous relaxation of the noisy discrete objective using only \(d+1\) function evaluations per simplex, providing gradient directions for free.

3. Convergence is almost sure to the set of \(N_1\)-local minima, and the probability of error decays exponentially in the sample size when the feasible region is finite.

4. The 10% sample size growth rule (\(m_{k+1} = \lceil 1.1 m_k \rceil\)) satisfies the convergence conditions easily. Most computation happens early when evaluations are cheap.

5. The KSL implementation follows the paper faithfully but adds practical limits (spliMaxIterations, lineSearchIterMax) and a confidence-interval-based stopping criterion not in the original paper.

6. R-SPLINE outperforms ISC and ISC-AHA in terms of oracle calls to near-optimality, especially as dimension grows, due to gradient-guided search and negligible overhead.

References

• Wang, H., Pasupathy, R., and Schmeiser, B. W. (2013). Integer-ordered simulation optimization using R-SPLINE: Retrospective Search with Piecewise-Linear Interpolation and Neighborhood Enumeration. ACM Transactions on Modeling and Computer Simulation, 23(3), Article 17.

• Hong, J. and Nelson, B. (2006). Discrete optimization via simulation using COMPASS. Operations Research, 54(1), 115–129.

• Xu, J., Hong, L. J., and Nelson, B. L. (2010). Industrial strength COMPASS: A comprehensive algorithm and software for optimization via simulation. ACM Trans. Model. Comput. Simul., 20, Article 19.

• Pasupathy, R. (2010). On choosing parameters in retrospective-approximation algorithms for stochastic root finding and simulation optimization. Operations Research, 58, 889–901.