Performative scenario optimization solutions converge to classical stochastic programming solutions as the strength of the decision-feedback effect decreases, with convergence rate bounded by the Lipschitz modulus of the distribution map.

This discovery has a Claude-generated validation package with a full experimental design.

Precise Hypothesis

Spine & Adversarial Read

• highThe claimed linear-in-L(θ) bound is asserted without citation to a specific theorem; performative prediction theory (Perdomo et al. 2020) establishes contraction/convergence of iterates under repeated risk minimization, but scenario optimization (sample-based, distributionally robust) has different finite-sample error interactions that may produce a different (e.g., sqrt or higher-order) rate — the EVP has not established this bound exists in prior literature, only that it is plausible by analogy.
  The methodology explicitly tests β (regression exponent) as an open empirical question rather than assuming β=1; success criteria allow β∈[0.8,1.2] as a falsifiable range. This is honest, but the EVP does not yet resolve whether the *scenario-optimization-specific* theoretical derivation exists — this must be attempted analytically (e.g., via McDiarmid/Wasserstein-based perturbation analysis of SAA solutions) before or alongside the empirical study, otherwise the empirical fit is descriptive, not confirmatory of a proven bound.
• mediumWhy these specific problem families (portfolio, newsvendor, strategic classification) and not others (e.g., queueing, energy dispatch)? The choice appears driven by convenience/availability of standard benchmarks rather than a principled argument that these span the relevant space of Lipschitz-modulus behaviors (linear vs. nonlinear feedback, convex vs. non-convex constraint sets).
  Partial justification given: these three families are the standard benchmarks in both performative prediction (strategic classification) and stochastic programming (portfolio, newsvendor) literatures, enabling direct comparison to prior published baselines. However, the EVP does not justify why 3 families and 3 perturbation classes constitute sufficient coverage of the hypothesis space — a more rigorous approach would include a fourth adversarially-designed 'stress-test' family specifically constructed to try to break the linear bound (e.g., highly non-smooth or discontinuous distribution maps), which is currently absent.
• mediumWasserstein distance approximation via Sinkhorn regularization introduces a bias that could systematically distort the measured L(θ), especially at small θ where the true signal is small — this could produce a spurious linear fit (regularization bias masquerading as true convergence rate) rather than reflecting genuine problem structure.
  Acknowledged in KNOWN_FAILURE_MODES and partially addressed by Checkpoint 3, but the EVP does not specify a concrete mitigation (e.g., exact OT solver cross-check via linear programming for small-to-moderate problem sizes, or entropic regularization parameter sweep to test sensitivity of β to Sinkhorn epsilon). This should be added as an explicit robustness check before final success/failure determination.

Experimental Protocol

ROI Projection

Implementation Sketch

# Pseudocode: Performative Scenario Optimization Convergence Test

def classical_SAA(P0, N=10000):
    scenarios = sample(P0, N)
    x0_star = solve_convex_program(scenarios)
    return x0_star

def performative_RRM(theta, P0, g, max_iter=50, tol=1e-6):
    x_t = initialize_x()
    for t in range(max_iter):
        P_theta_xt = P0 + theta * g(x_t)   # decision-dependent shift
        scenarios = sample(P_theta_xt, N=10000)
        x_next = solve_convex_program(scenarios)
        if norm(x_next - x_t) < tol:
            break
        x_t = x_next
    L_theta = wasserstein_distance(P_theta_xt, P0)  # Sinkhorn approx
    return x_t, L_theta

def run_experiment():
    x0_star = classical_SAA(P0)
    results = []
    for theta in [0, 0.01, 0.05, 0.1, 0.2, 0.4, 0.6, 0.8, 1.0]:
        for seed in range(30):
            x_theta, L_theta = performative_RRM(theta, P0, g, seed=seed)
            gap = norm(x_theta - x0_star)
            results.append((theta, L_theta, gap, seed))
    # regression: log(gap) ~ beta * log(L_theta) + log(C)
    beta, logC, R2 = fit_loglog(results)
    return results, beta, logC, R2