Repository navigation
[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
Description
Description
SyntheticControl solves for its donor weights, min_w ||x1 - X0 w||^2 over the simplex (with V folded in), by Frank-Wolfe. At the default inner_max_iter=10000, that solve often stops on its budget rather than its min-decrease rule on ordinary factor-model panels. In the reproducer's three designs, the placebo test failed closed or lost placebos on 40 of 40, 40 of 40 and 33 of 40 panels. Three things follow:
- the treated fit warns "did not converge", and
in_space_placebo()fails closed (placebo_p_valueNaN); - otherwise, placebos whose refit hit the budget are dropped from the reference set. Which units drop depends on how hard they are to fit, so the rank is taken against a selected subset rather than all units;
- in the default nested V search, evaluations that hit the budget are excluded from V ranking.
The budget flag does not separate good fits from bad ones. The inner problem has an exact finite solution (below). Measured against it, fits flagged at 10 000 steps were about as close to the optimum as unflagged ones. Twenty times the budget clears every flag, yet it barely moves the flagged fits and does not reach the optimum.
Minimal reproducer
V is held uniform (v_method="custom", standardize="none"), so only the inner solve is exercised. There are 40 panels per design.
"""diff-diff 3.12.0: SyntheticControl's Frank-Wolfe weight solve, at the default inner_max_iter,
exhausts its budget on a sizeable share of ordinary panels; the treated fit then fails closed
(placebo p-value NaN) or placebos are dropped as failed.
V is held uniform (v_method="custom", standardize="none") so only the inner simplex solve
min_w ||x1 - X0 w||^2 is exercised. Its exact optimum comes from one NNLS:
u = argmin_{u >= 0} || [X0 - x1 1'; 1'] u - [0; 1] ||, w = u / sum(u).
For u = s w the objective is s^2 a + (s - 1)^2 with a = ||X0 w - x1||^2; it is least at
s = 1 / (1 + a), where it equals a / (1 + a), increasing in a.
"""
import warnings
import numpy as np
import pandas as pd
from diff_diff import SyntheticControl
from scipy.optimize import nnls
def panel(rng, units, pre, post):
"""Unit level + three AR(1) factors with uniform loadings + AR(1) errors."""
periods = pre + post
factors = np.empty((periods, 3))
errors = np.empty((units, periods))
factors[0] = rng.normal(size=3) / np.sqrt(0.75)
errors[:, 0] = rng.normal(size=units) / np.sqrt(0.75)
for t in range(1, periods):
factors[t] = 0.5 * factors[t - 1] + rng.normal(size=3)
errors[:, t] = 0.5 * errors[:, t - 1] + rng.normal(size=units)
return rng.normal(size=(units, 1)) + rng.uniform(0, 1, (units, 3)) @ factors.T + errors
def long_frame(y, pre):
n, periods = y.shape
d = np.zeros((n, periods), dtype=int)
d[0, pre:] = 1
return pd.DataFrame({"unit": np.repeat(np.arange(n), periods), "time": np.tile(np.arange(periods), n),
"y": y.ravel(), "d": d.ravel()})
def exact_pre_rmspe(y, pre):
x1, x0 = y[0, :pre], y[1:, :pre].T
z = x0 - x1[:, None]
scale = np.sqrt(np.mean(z**2)) # keeps the two blocks comparable
u, _ = nnls(np.vstack([z / scale, np.ones(z.shape[1])]), np.r_[np.zeros(pre), 1.0], maxiter=5000)
w = u / u.sum()
return float(np.sqrt(np.mean((x1 - x0 @ w) ** 2)))
for inner_max_iter in (10_000, 200_000):
for units, pre, post in ((20, 20, 5), (40, 20, 5), (20, 40, 5)):
warned = nan_p = dropped = either = 0
excess = []
for seed in range(40):
y = panel(np.random.default_rng([units, pre, seed, 7]), units, pre, post)
with warnings.catch_warnings(record=True) as caught:
warnings.simplefilter("always")
fit = SyntheticControl(v_method="custom", custom_v=np.ones(pre), standardize="none",
inner_max_iter=inner_max_iter, seed=0).fit(
long_frame(y, pre), outcome="y", treatment="d", unit="unit", time="time",
pre_period_outcomes="all")
fit.in_space_placebo()
warned += any("did not converge" in str(w.message) for w in caught)
nan_p += not np.isfinite(fit.placebo_p_value)
dropped += int(fit.n_failed or 0) > 0
either += (not np.isfinite(fit.placebo_p_value)) or int(fit.n_failed or 0) > 0
excess.append(fit.pre_rmspe / exact_pre_rmspe(y, pre) - 1.0)
print(f"inner_max_iter={inner_max_iter}, {units - 1} donors, {pre} pre-periods, 40 panels: "
f"'did not converge' warning on {warned}, placebo_p_value NaN on {nan_p}, "
f"placebos dropped as failed on {dropped}, either on {either}; treated pre-RMSPE above the exact optimum "
f"by median {np.median(excess):.1e}, max {np.max(excess):.1e} (relative)", flush=True)Expected behavior
On ordinary panels, the weight solve should reach its optimum at the default settings, so that the placebo test keeps its full reference set.
Actual behavior
inner_max_iter=10000, 19 donors, 20 pre-periods, 40 panels: 'did not converge' warning on 10, placebo_p_value NaN on 10, placebos dropped as failed on 30, either on 40; treated pre-RMSPE above the exact optimum by median 1.2e-04, max 5.7e-04 (relative)
inner_max_iter=10000, 39 donors, 20 pre-periods, 40 panels: 'did not converge' warning on 15, placebo_p_value NaN on 15, placebos dropped as failed on 25, either on 40; treated pre-RMSPE above the exact optimum by median 1.9e-04, max 8.4e-04 (relative)
inner_max_iter=10000, 19 donors, 40 pre-periods, 40 panels: 'did not converge' warning on 3, placebo_p_value NaN on 3, placebos dropped as failed on 30, either on 33; treated pre-RMSPE above the exact optimum by median 1.2e-04, max 2.2e-04 (relative)
inner_max_iter=200000, 19 donors, 20 pre-periods, 40 panels: 'did not converge' warning on 0, placebo_p_value NaN on 0, placebos dropped as failed on 0, either on 0; treated pre-RMSPE above the exact optimum by median 1.1e-04, max 4.6e-04 (relative)
inner_max_iter=200000, 39 donors, 20 pre-periods, 40 panels: 'did not converge' warning on 0, placebo_p_value NaN on 0, placebos dropped as failed on 0, either on 0; treated pre-RMSPE above the exact optimum by median 1.6e-04, max 8.1e-04 (relative)
inner_max_iter=200000, 19 donors, 40 pre-periods, 40 panels: 'did not converge' warning on 0, placebo_p_value NaN on 0, placebos dropped as failed on 0, either on 0; treated pre-RMSPE above the exact optimum by median 1.2e-04, max 2.2e-04 (relative)Split by the flag, for the treated fit at 10 000 steps:
- 19 donors / 20 pre-periods: the flagged fits were a median 1.2e-4 above the exact optimum, the rest 1.2e-4. At 200 000 steps the flagged fits were 8.2e-5 above it.
- 39 donors: flagged 2.5e-4, the rest 1.6e-4; flagged at 200 000 steps, 1.7e-4.
- 19 donors / 40 pre-periods: flagged 1.4e-4, the rest 1.2e-4; flagged at 200 000 steps, 1.2e-4.
At the full defaults (v_method="nested", standardize="std"), the same inner solve runs at every V the search tries. Evaluations that hit the budget are excluded from V ranking, with the warning "Inner Frank-Wolfe did not converge on X of Y weight solves during nested V selection". I have not yet measured how often that happens on these panels, because each nested fit is slow. I will add the numbers here once the run finishes.
Proposed fix
With intercept=False, zeta=0, which is SyntheticControl's inner problem, the simplex-constrained least squares reduces exactly to one non-negative least squares:
Z = V^½ (X0 - x1 1') # columns: donor minus treated
u = argmin_{u >= 0} || [Z / s; 1'] u - [0; 1] ||^2 # s > 0 any scale, e.g. RMS of Z
w = u / sum(u)
The reduction works as follows. Write u = c w with w on the simplex and c = 1'u. The objective is c^2 a / s^2 + (c - 1)^2 with a = ||Z w||^2 = ||V^½ (X0 w - x1)||^2. Its minimum over c is a / (a + s^2), which increases with a. So the minimiser's w is the simplex least-squares solution, whatever s is. s only conditions the problem.
Lawson-Hanson's active set (scipy.optimize.nnls) terminates in finitely many steps. So there is no budget to tune, no convergence flag, and no residual gap. It is what the reproducer uses as exact_pre_rmspe.
SyntheticDiD's regularised and intercepted Frank-Wolfe, which matches R's synthdid, would be untouched. If Frank-Wolfe is preferred for SyntheticControl too, away-step or pairwise Frank-Wolfe (Lacoste-Julien & Jaggi 2015) converge linearly over a polytope. Classic Frank-Wolfe converges sublinearly when the optimum lies on a face of the simplex, which is the usual case for synthetic control weights.
Version
diff-diff 3.12.0 (PyPI wheel, Rust backend), Python 3.14.7, Linux.
Activity
ilgrad commentedon Sep 30, 2026
The numbers at the full defaults (
v_method="nested",standardize="std"), for the treated fit only, on 10 panels per design from the reproducer'spanel():None of the 30 final fits warned "did not converge". So at the defaults, the budget acts mostly inside the V search. There, V is ranked only among the candidates whose inner solve stopped on its min-decrease rule, and the only sign of that is the aggregate warning.
With the NNLS reduction, every inner solve returns the exact optimum, so no evaluation would be excluded.
ilgrad commentedon Oct 5, 2026
A patch for the reduction proposed above, what it changes, and the choices it leaves to you.
The patch, against
mainat 6874b5d.nnlsis SciPy's, already a dependency; it raises when it runs out of iterations, so a solve that does not finish fails loudly rather than returning a truncatedW.diff
On the reproducer's panels (uniform V, 40 panels in each of the three designs): no fit warns that it did not converge,
placebo_p_valueis finite on every panel, and no placebo is dropped. The treated fits' pre-RMSPE equals the exact optimum to 4e-16, relative. As an independent check, the simplex QP's KKT conditions hold at every returnedW: stationarity on the support within 1.3e-15 of the largest gradient entry, and no dual infeasibility off it.On the test suite (
tests/test_methodology_synthetic_control.pyandtests/test_t25_synthetic_control_policy_drift.py,DIFF_DIFF_BACKEND=python): 277 passed onmain; with the patch 273 passed and 4 failed.test_inner_v_search_nonconvergence_warningandtest_cv_fails_closed_when_training_solve_truncates[specs0]and[specs1]force a truncated inner solve withinner_max_iter=1. The solve no longer truncates.test_top_donor_weights(Tutorial 25): unit 3's weight is 0.414, outside [0.20, 0.32].Why Tutorial 25's weights move. With three predictors and 20 donors, the treated unit's predictors lie inside the donors' convex hull. The inner problem's optimum is then zero, attained on a whole face of the simplex. At the V
mainselects, Frank-Wolfe (at the tutorial'sinner_min_decrease=1e-3) stops at an inner loss of 2.7e-4, and the exact solve reaches 1e-33, with weights up to 0.157 apart. So the nested search ranks each V by which zero-lossWthe solver happens to return. With the patch the tutorial's fit has a pre-RMSPE of 0.934 against 0.987 onmain, an ATT of 6.69 against 6.79, and weights of 0.497 and 0.414 on units 2 and 3 against 0.509 and 0.257.What the patch leaves to you:
W. NNLS returns a vertex of the optimal face: sparse and deterministic. Frank-Wolfe returns wherever its budget left it, and an interior-point solver such as RSynth'sipoptends to a point inside the face. If the reported weights should follow a stated rule there, a tie-break can be added: Abadie and L'Hour's (2021) penalty, or the minimum-normWamong the optima.inner_max_iterandinner_min_decreaseno longer act on SyntheticControl's weights. They can be deprecated, or kept, accepted and ignored, with a note.in_space_placebo()failing closed or dropping placebos, and the CV fail-closed. Their tests can keep their coverage by monkeypatching_inner_solve_Wto report a truncated solve, as one test in the file already does.conformal.py:108(CWZ's proxy fit,zeta=0) andprep.py:1053solve the same simplex least squares by Frank-Wolfe. Withzeta > 0the ridge term stacks a multiple of the identity under the predictors, and the same reduction applies. The patch leaves both alone.The alternative that keeps Frank-Wolfe is its away-step or pairwise variant (Lacoste-Julien and Jaggi 2015), which converges linearly on a polytope. It keeps the parameters and the non-convergence paths meaningful. It needs both the NumPy and the Rust kernels changed, still stops on a budget, and faces the same question of ties.
I can open the PR for either: the NNLS solve, with 1–3 as you choose, or away-step Frank-Wolfe.
ilgrad commentedon Oct 6, 2026
PR #842 implements the NNLS route from my comment above. Two things changed from that patch:
scipy.optimize.nnls. At an exact-fit tie the vertex an active set returns depends on its pivoting, and SciPy's implementation has changed across releases: on Tutorial 25, SciPy 1.11 to 1.18 gave five different sets of placebo RMSPE ratios, and 1.15 an unbounded constant-effect set. With the in-house solve they agree to the sixth decimal on 1.11 to 1.18.Vto a corner on Tutorial 25),inner_max_iter/inner_min_decreaseaccepted but unused (a TODO row), the unreachable non-convergence paths kept and tested through a monkeypatched truncation, and both siblings (the CWZ proxy,rank_control_units'synthetic_weight) solved the same way.The PR body has the fail-before run, the SciPy 1.10 to 1.18 matrix, and every number that moves.