Skip to content

fix(synthetic_control): solve the donor weights exactly by non-negative least squares - #842

Open
ilgrad wants to merge 1 commit into
igerber:mainfrom
ilgrad:fix/sc-exact-simplex-weights
Open

ilgrad wants to merge 1 commit into
igerber:mainfrom
ilgrad:fix/sc-exact-simplex-weights

Conversation

@ilgrad

@ilgrad ilgrad commented Oct 6, 2026

Copy link
Copy Markdown

Fixes #838.

Problem

SyntheticControl solved its donor weights, min_W ||V^½ (X1 − X0 W)||² over the simplex, by Frank-Wolfe with a budget of inner_max_iter=10000 steps. Frank-Wolfe converges sublinearly when the optimum lies on a face of the simplex, the usual case for synthetic control weights, and on ordinary factor-model panels it stopped on that budget. The treated fit then warned "did not converge" and the diagnostics that need a converged fit were skipped (placebo_p_value NaN); placebos whose refit hit the budget were dropped from the reference set; and the nested V search ranked V only among the solves that met their stopping rule. In #838's reproducer the placebo test failed closed or lost placebos on 40, 40 and 33 of 40 panels, and at the full defaults the nested search excluded 46%, 25% and 0.5% of its solves from V ranking.

Fix

The problem reduces exactly to one non-negative least squares (the derivation is in #838 and in _sc_weight_nnls's docstring): with Z = V^½ (X0 − X1 1'), u = argmin_{u ≥ 0} ||[Z / s; 1'] u − [0; 1]|| and W = u / 1'u, for any scale s > 0. An active-set method solves it in finitely many steps, so there is no budget to exhaust.

  • utils._sc_weight_nnls(Y, zeta): _sc_weight_fw's uncentered problem (intercept=False), zeta² ||w||² + ||A w − b||² / N on the simplex, solved exactly. zeta > 0 stacks √N ζ I under Z.
  • utils._nnls_lawson_hanson(M, d, max_iter): Lawson and Hanson's active set (1974, ch. 23), with their guard against an entering column that rounding leaves without positive weight. It is ~40 lines of NumPy rather than a call to scipy.optimize.nnls; see SciPy versions. Past max_iter (50 per donor) it raises RuntimeError rather than return a truncated solve; no handler in these modules catches it, and no test or run below raised it.
  • _inner_solve_W, which every V evaluation and every placebo, leave-one-out, in-time and sparse refit goes through, calls _sc_weight_nnls and always reports converged.
  • Two siblings solved the same problem by Frank-Wolfe and now use the same solve: CWZ's conformal proxy (conformal._cwz_proxy_fit, zeta=0), and prep.rank_control_units's informational synthetic_weight column (zeta = sqrt(lambda_reg / N), which stopped at max_iter=1000). A panel with a missing pre-period outcome keeps the uniform column it had.
  • Warnings and docstrings that named Frank-Wolfe, or prescribed inner_max_iter / inner_min_decrease as the remedy, now name the outer search's knobs (optimizer_options['maxiter'], n_starts), or none where only the inner solve could fail.
  • SyntheticDiD's Frank-Wolfe (intercept, ridge, R synthdid parity) is untouched.

Choices

The four from my comment on #838:

  1. Ties. Where the treated predictors lie inside the donors' convex hull, every W on a face of the simplex fits them exactly. The active set returns a vertex of that face (at most k + 1 donors); R Synth's ipop tends to a point inside it, and Frank-Wolfe stopped wherever its budget left it. REGISTRY records this as a **Note (deviation from R)**. I added no selection rule. I tried the minimum-norm W (a 1e-10 relative ridge) on Tutorial 25: W(V) becomes continuous, the nested search then drives V to (0.017, 2e-12, 0.983), a corner where one predictor's weight is negligible, and the ATT to 6.91. A TODO row records the question; Abadie and L'Hour (2021) is the published rule.
  2. inner_max_iter, inner_min_decrease are accepted and validated but no longer used, and their docstrings say so. Deprecating them is a TODO row.
  3. The inner non-convergence paths stay but cannot fire for the inner solve: the nested and CV aggregate warnings, the fail-closed and dropped-placebo branches, and the conformal indeterminate grid points. Their tests keep their coverage by monkeypatching _inner_solve_W to report a truncated solve (_truncate_inner_solves). Keeping or removing the paths is in the same TODO row.
  4. Siblings: both solved exactly here (above).

SciPy versions

At a tie, the vertex an active set returns depends on its pivoting, and scipy.optimize.nnls has been reimplemented several times. With it, Tutorial 25's 21 placebo RMSPE ratios took five different sets over SciPy 1.11 to 1.18, where main's Frank-Wolfe ratios agree to 2e-5, and under 1.15 the constant-effect confidence set was unbounded. With _nnls_lawson_hanson the ratios agree to the sixth decimal under 1.11 to 1.18.

Tutorial 25 under each release, DIFF_DIFF_BACKEND=python, Python 3.11 (1.17.1 on 3.12, 1.18.1 on 3.14). A letter names one set of the 21 placebo RMSPE ratios, compared at six decimals; the brackets are the constant-effect set as the tutorial prints it.

SciPy NumPy main this PR this PR with scipy.optimize.nnls
1.10.1 1.26.4 A′ [4.16, 24.44] B′ [3.53, 31.50] C1 [3.39, 31.50]
1.11.4 1.26.4 A [4.16, 24.44] B [3.53, 31.50] C2 [3.39, 31.50]
1.12.0 1.26.4 A [4.16, 24.44] B [3.53, 31.50] C3 [3.39, 31.36]
1.13.1 2.0.2 A [4.16, 24.44] B [3.53, 31.50] C3 [3.39, 31.36]
1.14.1 2.1.3 A [4.16, 24.44] B [3.53, 31.50] C3 [3.39, 31.36]
1.15.3 2.2.6 A [4.16, 24.44] B [3.53, 31.50] C4, unbounded
1.16.3 2.3.5 A [4.16, 24.44] B [3.53, 31.50] C5 [3.53, 31.50]
1.17.1 2.4.6 A* [4.05, 24.44] B [3.39, 31.50] C6 [3.53, 31.50]
1.18.1 2.5.3 A* [4.05, 24.44] B [3.39, 31.36] C6 [3.53, 31.50]

A* is A with one ratio 1.5e-5 away (1.119307 → 1.119322). A′ and B′ each move one placebo under SciPy 1.10; that comes from the outer V search, not the inner solve, on main as on this branch.

The constant set's printed endpoints move by one grid step between environments on both branches (main: 4.16 in the first seven rows, 4.05 in the last two). That is older than this PR: the tutorial prints the first and last accepted points of the 200-point display grid that _invert_sharp_null lays between the exact endpoints, and at those two points a placebo ties the treated unit, so rounding decides whether they show as accepted. The exact endpoints in effect_confidence_set do not depend on that tie. I left it alone here.

Behavioral changes

  • Estimates move by the former solver's gap. On [Bug]: SyntheticControl's Frank-Wolfe stops on its budget at the default inner_max_iter; in_space_placebo() then fails closed or drops placebos on 33-40 of 40 panels #838's panels the treated pre-RMSPE was a median 1e-4 above the optimum (relative); it now equals it to 7e-16. At a tie the move can be large, because Frank-Wolfe's stopping point and the vertex are different exact fits. Tutorial 25 is such a case, and its text and drift bands are updated:

    Tutorial 25 main this PR
    weights 2: 0.509, 3: 0.257, 13: 0.082, 18: 0.071, 16: 0.016, and 0.0043 on each of 15 more 2: 0.497, 3: 0.414, 4: 0.076, 8: 0.014
    largest predictor gap, synthetic vs treated 0.037 0.000
    ATT 6.79 6.69
    pre-RMSPE 0.987 0.934
    placebo p 1/21 1/21
    constant-effect set (as printed) [4.05, 24.44] [3.39, 31.36]
    leave-one-out, max |ΔATT| 0.76 (20 refits) 0.50 (4 refits)
    in-time placebo ATTs (backdated 20, 35, 50) 0.598, 0.075, −0.367 0.150, 0.547, 0.916
    CWZ joint p, true path 0.831 0.708
    CWZ joint p, zero 0.015 0.015
    average-effect CI [6.50, 7.48] [6.46, 7.79]

    The tutorial's prose still holds (every in-time placebo below 1 in absolute value, the leave-one-out bound, the conformal claims); the numbers quoted in it are updated.

  • Conformal intervals are narrower and still cover above nominal. Over 50 simulated panels (6 donors, 22 periods, constant effect 6; 300 period intervals at nominal 90%, inverted on a fine grid), coverage goes from 0.970 to 0.953 and the median width from 0.515 to 0.450. No grid point failed to converge on either branch.

  • One existing test now inverts on a fine grid: test_conformal_confidence_intervals_recover_true_constant_effect. The default grid puts 100 points on ±0.5·|estimate|, a step of about 0.06 there, against intervals about 0.2 wide, so each interval's endpoints sit up to a step inside the accepted set. On that grid this branch brackets the true effect in 4 of 6 periods (the test asks for 5); on bounds=(5, 7), n_grid=401 it brackets 5 of 6, and main 6 of 6 on both grids. A TODO row records the grid floor.

  • R parity. At R's V (Basque Tier 1), the weights match Synth to 4.5e-6 (9.6e-5 on main); the test's tolerance goes from 1e-3 to 2e-5. At the defaults (nested V, seed=0), the pre-RMSPE is 0.0888 against Synth's 0.0942 and the ATT −0.5842 against −0.5799, with weights 0.8331 (Cataluña) and 0.1669 (Madrid).

Methodology references (required if estimator / math changes)

  • Method name(s): the synthetic control weight problem (Abadie, Diamond & Hainmueller 2010), CWZ's conformal proxy (Chernozhukov, Wüthrich & Zhu 2021); non-negative least squares by Lawson and Hanson's active set (Solving Least Squares Problems, 1974, ch. 23).
  • Paper / source link(s): https://doi.org/10.1198/jasa.2009.ap08746, https://doi.org/10.1080/01621459.2021.1920957
  • Any intentional deviations from the source (and why): none from the papers; the estimator is the same, solved to its optimum. From R Synth: at an exact-fit tie the reported W is the active set's vertex where ipop tends to an interior point (REGISTRY Note, above).

Validation

  • Tests added/updated.
    • New, test_methodology_synthetic_control.py: _sc_weight_nnls meets the simplex QP's KKT conditions (3 factor panels, zeta 0 and 0.3) and returns a vertex at an exact-fit tie; _nnls_lawson_hanson equals scipy.optimize.nnls where the optimum is unique, meets KKT on tall, square and wide matrices with a duplicated column, skips a column that rounding lets in, and raises rather than truncate; on factor panels where the budgeted solve stopped, the inner solve is exact, in_space_placebo() keeps every placebo, and the conformal proxy converges.
    • New, test_prep.py: synthetic_weight is the exact simplex optimum, and stays uniform on a panel with a missing pre-period outcome.
    • Updated: the six non-convergence path tests now truncate through _truncate_inner_solves (the solve itself cannot fail); the CWZ proxy oracle is held to SLSQP at atol=1e-7 with a one-step budget (was 1e-2 with 200 000 steps); Basque Tier 1 at atol=2e-5 (was 1e-3); the conformal-interval test on a fine grid (below); Tutorial 25's drift bands; the naming guard's zeta allowlist.
  • Fail before, pass after. The branch's three test files against main's code (6874b5d), same environment: 32 failed, 526 passed.
    • 23 are the new solver tests: _sc_weight_nnls and _nnls_lawson_hanson do not exist on main.
    • 7 fail on main's results: test_inner_solve_is_exact_where_a_budgeted_solve_stopped (the budgeted solve does not converge); test_in_space_placebo_keeps_every_placebo_on_a_factor_panel[0] and [1] (the treated fit warns "did not converge", and 6 of 19 placebos are dropped); test_conformal_proxy_converges_on_a_factor_panel; test_cwz_proxy_fit_matches_scipy_simplex_ls; test_synthetic_weight_is_the_exact_simplex_optimum (main's synthetic_weight objective 8.9e-11 against the optimum's 3.1e-14); and test_basque_tier1_custom_v_parity at atol=2e-5.
    • 2 fail on what this PR changes on purpose: test_top_donor_weights (Tutorial 25's weights) and test_inner_nonconvergence_warning (the warning no longer names Frank-Wolfe).
  • [Bug]: SyntheticControl's Frank-Wolfe stops on its budget at the default inner_max_iter; in_space_placebo() then fails closed or drops placebos on 33-40 of 40 panels #838's reproducer on this branch, both budgets (inner_max_iter no longer acts): no "did not converge" warning, no NaN p-value and no dropped placebo on any of the 120 panels, and the treated pre-RMSPE above the exact optimum by at most 6.7e-16 (relative). The nested count from my comment: no solve excluded from V ranking on any of the 30 fits.
  • Commands: Python 3.14.7, numpy 2.5.3, pandas 3.0.6, scipy 1.18.1, DIFF_DIFF_BACKEND=python.
    • pytest -n 4 -m '' (slow tests included) over test_methodology_synthetic_control.py, test_t25_synthetic_control_policy_drift.py, test_prep.py, test_utils.py, test_methodology_sdid.py, test_diagnostic_report.py, test_business_report.py, test_practitioner.py, test_rust_backend.py, test_aliases.py, test_v4_matrix.py, test_v4_wrapper_shims.py, test_naming_guard.py, test_changelog_fragments.py, test_guides.py, test_doc_snippets.py, test_doc_deps_integrity.py, test_docs_ia.py and test_tracking_files.py: 2324 passed, 221 skipped, 2 failed.
      • test_naming_guard.py::test_allowlists_are_reachable: this PR removed the last zeta= from synthetic_control.py, which left its zeta allowlist entry with nothing to cover. The entry is removed; the docs, changelog, prep, guides and naming-guard files then pass (899 passed, 7 skipped).
      • test_rust_backend.py::TestClusterVcovDeterminism::test_raw_kernel_noncontiguous_ids_bit_identical fails on main too here: that class has no skipif(not HAS_RUST_BACKEND) and imports diff_diff._rust_backend, which is not built here.
    • ruff check, black --check and mypy on the changed modules: clean. changelog_compile.py check: OK.
    • Tutorial 25 re-executed: no errors; its drift tests pass.
  • SciPy floor: the solver tests, the T25 drift tests and the synthetic_weight tests (50, selected with -k "t25 or nnls or lawson or exact or synthetic_weight or factor_panel" over the three files) pass under SciPy 1.10.1, 1.11.4, 1.12.0, 1.13.1, 1.14.1, 1.15.3, 1.16.3 (Python 3.11) and 1.17.1 (Python 3.12), NumPy as in the table above.
  • Not run: the Rust backend (no toolchain here; its tests skip but for the one above), the rest of the suite, and the other notebooks.

Security / privacy

  • Confirm no secrets/PII in this PR: confirmed. Code, tests, docs and one re-executed tutorial; the simulations use synthetic data.

Changelog

  • changelog.d/ fragment added (or N/A - no user-visible change): changelog.d/20261006-sc-exact-simplex-weights.md, ### Fixed and ### Behavioral Changes. It states which results move and by how much (Tutorial 25, Basque, the conformal coverage), the vertex at ties, and that inner_max_iter / inner_min_decrease no longer act.

…ve least squares

Frank-Wolfe with a budget of inner_max_iter steps stopped short of the simplex
least-squares optimum on ordinary factor-model panels, so the treated fit failed
closed, placebos were dropped as failed, and the nested V search ranked V by where
truncated solves stopped. The problem reduces exactly to one non-negative least
squares, solved here by Lawson-Hanson's active set, which terminates finitely and,
being diff-diff's own, returns the same vertex at a tie under every SciPy release.
The CWZ conformal proxy and rank_control_units' synthetic_weight column solve the
same problem and use the same solve.

Fixes igerber#838

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

1 participant