[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

[ilgrad]

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_value NaN);
• 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.

[ilgrad]

The numbers at the full defaults (v_method="nested", standardize="std"), for the treated fit only, on 10 panels per design from the reproducer's panel():

donors	pre-periods	fits whose V search excluded solves	solves excluded from V ranking
19	20	7 of 10	81 151 of 176 080 (46%)
39	20	7 of 10	31 968 of 128 176 (25%)
19	40	3 of 10	1 097 of 216 987 (0.5%)

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.

import re
import warnings

import numpy as np
from diff_diff import SyntheticControl

# panel() and long_frame() as in the reproducer above
for units, pre in ((20, 20), (40, 20), (20, 40)):
    fits = excluded = solves = final = 0
    for seed in range(10):
        y = panel(np.random.default_rng([units, pre, seed, 7]), units, pre, 5)
        with warnings.catch_warnings(record=True) as caught:
            warnings.simplefilter("always")
            SyntheticControl(seed=0).fit(long_frame(y, pre), outcome="y", treatment="d", unit="unit",
                                         time="time", pre_period_outcomes="all")
        for w in caught:
            m = re.search(r"Inner Frank-Wolfe did not converge on (\d+) of (\d+)", str(w.message))
            if m:
                fits += 1
                excluded += int(m.group(1))
                solves += int(m.group(2))
            elif "did not converge" in str(w.message):
                final += 1
    print(units - 1, pre, fits, excluded, solves, final)

With the NNLS reduction, every inner solve returns the exact optimum, so no evaluation would be excluded.

[ilgrad]

A patch for the reduction proposed above, what it changes, and the choices it leaves to you.

The patch, against main at 6874b5d. nnls is 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 truncated W.

diff

diff --git a/diff_diff/synthetic_control.py b/diff_diff/synthetic_control.py
index bc61afb..ac64e70 100644
--- a/diff_diff/synthetic_control.py
+++ b/diff_diff/synthetic_control.py
@@ -38,14 +38,14 @@ from typing import Any, Dict, List, Optional, Tuple, cast
 
 import numpy as np
 import pandas as pd
-from scipy.optimize import minimize
+from scipy.optimize import minimize, nnls
 
 from diff_diff._base import BaseEstimator
 from diff_diff.synthetic_control_results import (
     SyntheticControlResults,
     _SyntheticControlFitSnapshot,
 )
-from diff_diff.utils import _sc_weight_fw, safe_inference, warn_if_not_converged
+from diff_diff.utils import safe_inference, warn_if_not_converged
 
 __all__ = ["SyntheticControl", "synthetic_control"]
 
@@ -1141,7 +1141,12 @@ def _inner_solve_W(
 ) -> Tuple[np.ndarray, bool]:
     """Solve W*(V) = argmin_W (X1-X0 W)' diag(V) (X1-X0 W) on the unit simplex.
 
-    Folds V^½ into the predictors and delegates to the Frank-Wolfe simplex solver.
+    Solved exactly by one non-negative least squares: with ``Z = V^½ (X0 - X1 1')``,
+    ``u = argmin_{u >= 0} ||[Z / s; 1'] u - [0; 1]||`` and ``W = u / 1'u``. For
+    ``u = c W`` the objective is ``c² ||Z W||² / s² + (c - 1)²``, whose least value
+    over ``c`` increases with ``||Z W||²``, so ``W`` is the simplex optimum for any
+    scale ``s > 0``. Lawson-Hanson's active set terminates finitely, so the solve
+    always converges and ``inner_max_iter`` / ``inner_min_decrease`` do not apply.
     Returns the RAW simplex weights (no sparsification) so the outer objective is
     not perturbed; reporting-level cleanup happens in ``fit``.
     """
@@ -1149,21 +1154,13 @@ def _inner_solve_W(
     if J == 1:
         return np.array([1.0]), True
     vh = np.sqrt(v)
-    A = vh[:, None] * X0s  # (k, J)
-    b = vh * X1s  # (k,)
-    packed = np.column_stack([A, b])  # (k, J+1); last column is the target
-    min_decrease = inner_min_decrease * max(float(np.linalg.norm(b)), 1e-12)
-    with warnings.catch_warnings():
-        warnings.filterwarnings("ignore", message=r".*did not converge.*", category=UserWarning)
-        w, converged = _sc_weight_fw(
-            packed,
-            zeta=0.0,
-            intercept=False,
-            min_decrease=min_decrease,
-            max_iter=inner_max_iter,
-            return_convergence=True,
-        )
-    return np.asarray(w, dtype=float), bool(converged)
+    Z = vh[:, None] * (X0s - X1s[:, None])  # (k, J): donor minus treated
+    scale = float(np.sqrt(np.mean(Z**2))) or 1.0
+    A = np.vstack([Z / scale, np.ones(J)])
+    b = np.zeros(A.shape[0])
+    b[-1] = 1.0
+    u, _ = nnls(A, b, maxiter=50 * J)
+    return u / u.sum(), True
 
 
 def _v_starts(

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_value is 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 returned W: 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.py and tests/test_t25_synthetic_control_policy_drift.py, DIFF_DIFF_BACKEND=python): 277 passed on main; with the patch 273 passed and 4 failed.

• test_inner_v_search_nonconvergence_warning and test_cv_fails_closed_when_training_solve_truncates[specs0] and [specs1] force a truncated inner solve with inner_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 main selects, Frank-Wolfe (at the tutorial's inner_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-loss W the solver happens to return. With the patch the tutorial's fit has a pre-RMSPE of 0.934 against 0.987 on main, 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:

1. Ties in 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 R Synth's ipop tends 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-norm W among the optima.
2. inner_max_iter and inner_min_decrease no longer act on SyntheticControl's weights. They can be deprecated, or kept, accepted and ignored, with a note.
3. The non-convergence paths become unreachable for SyntheticControl: the nested search's aggregate warning, in_space_placebo() failing closed or dropping placebos, and the CV fail-closed. Their tests can keep their coverage by monkeypatching _inner_solve_W to report a truncated solve, as one test in the file already does.
4. Siblings. conformal.py:108 (CWZ's proxy fit, zeta=0) and prep.py:1053 solve the same simplex least squares by Frank-Wolfe. With zeta > 0 the 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]

PR #842 implements the NNLS route from my comment above. Two things changed from that patch:

• The NNLS is an in-house Lawson-Hanson active set (~40 lines of NumPy), not 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.
• The four choices are as I proposed them: no tie rule (a TODO row; a minimum-norm tie-break drives the nested V to a corner on Tutorial 25), inner_max_iter / inner_min_decrease accepted 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.