[Mirror] SC exact NNLS donor weights - CI-review mirror of PR #842

TODO.md

@@ -64,6 +64,8 @@ Related tracking surfaces:
 | Recover the Supplementary Appendix DGPs 2 / 3 / 5 (nonlinear `W` / `X` terms) on the bad-control ML lane - needs a flexible-learner fixture via the duck-typed learner route (the parametric built-ins are inconsistent there by design, Table S1). | `tests/test_methodology_dml_did.py` | bad-controls PR-B | Mid | Low |
 | DurationDiD alternative finite-sample PH estimators (the repaired cumulative-increment slope of Equation 3.5 and the average-hazard slope of Equation 3.6) as an opt-in beside the shipped mean-of-ratios coefficient; they coincide under exact PH and differ under sampling noise (review table, lines 306-311) | `duration_did.py` | #822 | Mid | Low |
 | `DiagnosticReport` / `BusinessReport` admission of `DurationDiDResults`: both reject it by type today (their batteries are keyed to mean-outcome parallel-trends diagnostics); admission needs a hazard-restriction check family (the Algorithm 2 pretest, the invalid-curve statuses, per-family bootstrap failure counts) and a narrative for per-date absorption ATTs | `diagnostic_report.py`, `business_report.py` | #822 | Mid | Low |
+| Conformal `_auto_grid`'s half-width floor `0.5 * abs(center)` dominates `8 * spread` when the effect is large relative to the residual scale: a narrow interval then holds only 3-4 of the 100 default grid points, and its endpoints move by a grid step. Size the grid from the residual scale, or refine around the accepted set. | `conformal.py` | #838 | Quick | Low |
+| `SyntheticControl` at an exact-fit tie (treated predictors inside the donors' convex hull): `W` and `V` are not identified, and the inner solve reports the active set's vertex (REGISTRY §SyntheticControl). A stated selection rule (Abadie-L'Hour 2021 penalty, or the minimum-norm `W`) is a methodology choice; note that any strictly convex tie-break makes `W(V)` continuous, and the nested search then moves `V` toward a corner where a predictor's weight is negligible (Tutorial 25 under a 1e-10 relative ridge: `V` = (0.017, 2e-12, 0.983)). | `synthetic_control.py`, `utils.py` | #838 | Mid | Low |
 
 ### Performance
 

changelog.d/20261006-sc-exact-simplex-weights.md

@@ -0,0 +1,56 @@
+### Fixed
+- **`SyntheticControl` donor weights are the exact simplex least-squares
+  solution** (#838): the inner weight solve, run at every `V` the nested
+  search tries and in every placebo refit, was Frank-Wolfe with an iteration
+  budget (`inner_max_iter=10000`) that ordinary factor-model panels often
+  exhausted. 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 search ranked `V` only among the solves that met their stopping
+  rule. The problem reduces exactly to one non-negative least squares, solved
+  by Lawson-Hanson's active set, which terminates in finitely many steps, so
+  the solve now reaches the optimum and none of those paths fires for it (a
+  solve that exhausted the active set's step guard, not reached in practice,
+  would return a feasible iterate flagged non-converged and take the same
+  fail-closed paths as before). On the issue's
+  three designs, 40 panels each, the placebo test had failed closed or lost
+  placebos on 40, 40 and 33 panels; it now keeps every placebo on all 120, and
+  the treated pre-RMSPE equals the exact optimum to rounding.
+  - The same solve replaces Frank-Wolfe in the conformal proxy
+    (`conformal_test()`, `conformal_average_effect()` and
+    `conformal_confidence_intervals()`: `proxy_converged` is true and
+    `n_grid_nonconverged` is 0 in practice) and in the informational `synthetic_weight`
+    column of `rank_control_units()`, which stopped at 1000 iterations; that
+    column's values move, the ranking does not. A panel with a missing
+    pre-period outcome keeps its uniform `synthetic_weight` column; a solve that
+    does not converge (not reached in practice) leaves the column NaN with a
+    `UserWarning` instead of presenting a fallback as the solution.
+
+### Behavioral Changes
+- **`SyntheticControl` results move by the former solver's gap** (#838): the
+  donor weights, `att`, `pre_rmspe`, and the placebo, leave-one-out, in-time,
+  sparse and conformal results. On generic panels the change is small: the
+  issue's treated pre-RMSPE was a median 1e-4 above the optimum, relative.
+  Where the treated unit's predictors lie inside the donors' convex hull,
+  every `V` fits them exactly, so the predictors identify neither the weights
+  nor `V`. The solve then returns a vertex of the optimal face; Frank-Wolfe
+  stopped wherever its budget left it, and R `Synth`'s interior-point solver
+  tends to a point inside the face. The vertex is the same under every SciPy
+  version: the active set is diff-diff's own, because `scipy.optimize.nnls`
+  returns a different vertex at a tie from one release to another (five
+  different sets of T25's placebo ratios over SciPy 1.11 to 1.18). The T25
+  tutorial panel is such a case: its weights go from five donors plus 0.0043
+  on each of the other 15 to four donors, `att` from 6.79 to 6.69 and
+  `pre_rmspe` from 0.987 to 0.934. On the Basque panel at R's `V`, the weights
+  match R `Synth` to 4.5e-6 (9.6e-5 before).
+- **Conformal p-values and intervals move** (#838): with the exact proxy the
+  intervals are narrower on simulated panels and still cover above nominal.
+  Over 50 panels (300 period intervals at nominal 90%), coverage goes from
+  0.970 to 0.953 and the median width from 0.515 to 0.450.
+
+### Deprecated
+- **`SyntheticControl(inner_max_iter=, inner_min_decrease=)`** (#838, ledger
+  M-148 / M-149): both now default to `None` and emit a `DeprecationWarning` when
+  passed; the value is discarded (`get_params()` reports `None`). The inner
+  weight solves have no iteration budget or stopping tolerance to set. Drop the
+  arguments; they are removed in 5.0.

diff_diff/conformal.py

@@ -26,17 +26,16 @@
 same objective when ``covariates=`` is supplied (the paper's note after eq 6),
 with outcome-only residuals. This is DISTINCT from the headline ADH V-matrix weights — CWZ's
 exactness theory (Lemma 1, Appendix D exchangeability) requires a time-symmetric
-proxy, which the ADH pre-period V-fit is not. Reuses the Frank-Wolfe simplex
-solver :func:`diff_diff.utils._sc_weight_fw`.
+proxy, which the ADH pre-period V-fit is not. The proxy is solved exactly, by one
+non-negative least squares (:func:`diff_diff.utils._sc_weight_nnls`).
 """
 
 import math
-import warnings
 from typing import Any, Dict, List, Optional, Tuple
 
 import numpy as np
 
-from diff_diff.utils import _sc_weight_fw
+from diff_diff.utils import _sc_weight_nnls
 
 __all__: List[str] = []
 
@@ -52,24 +51,23 @@ def _cwz_proxy_fit(
     y1: np.ndarray,
     Y0: np.ndarray,
     *,
-    max_iter: int,
-    min_decrease: float,
-    init_weights: Optional[np.ndarray] = None,
     x1_rows: Optional[np.ndarray] = None,
     X0_rows: Optional[np.ndarray] = None,
 ) -> Tuple[np.ndarray, np.ndarray, bool]:
     """Fit the canonical CWZ constrained-LS SC proxy over ALL given periods (eq 4).
 
     Minimises ``Σ_t (y1_t − Σ_j w_j Y0_{tj})²`` s.t. ``w ≥ 0, Σ w = 1`` (no
     intercept, no V-matrix). ``y1`` is the (already null-imputed) treated outcome
-    vector ``(T,)``; ``Y0`` the donor outcomes ``(T, J)``. Reuses
-    :func:`_sc_weight_fw` packed ``(T, J+1)`` (donor columns first, target last),
-    ``zeta=0, intercept=False`` — the solver projects onto the unit simplex, so the
-    ``w ≥ 0, Σ w = 1`` constraint is delivered with no extra normalization.
+    vector ``(T,)``; ``Y0`` the donor outcomes ``(T, J)``. Solved exactly by
+    :func:`_sc_weight_nnls` on the packed ``(T, J+1)`` (donor columns first, target
+    last), ``zeta=0`` — one non-negative least squares whose solution lies on the unit
+    simplex, so the ``w ≥ 0, Σ w = 1`` constraint is delivered with no extra
+    normalization.
 
-    ``min_decrease`` is the ALREADY-SCALED absolute convergence tolerance; the
-    caller must pass a θ0-invariant scale (e.g. the pre-window outcome norm) so CI
-    membership does not drift with the grid value via the tolerance.
+    The solve has no iteration budget to tune; ``converged`` is False only if the
+    active set exhausted its step guard or the data left the float64 range (neither
+    reached in practice), in which case the
+    caller treats the grid point as indeterminate.
 
     **Covariates (CWZ §2.3.2, after eq 6):** "it is straightforward to
     incorporate (transformations of) covariates X_jt into the estimation
@@ -103,17 +101,7 @@ def _cwz_proxy_fit(
         w = np.array([1.0], dtype=float)
         return w, y1 - Y0 @ w, True
     packed = np.column_stack([Y0_fit, y1_fit])  # (T+R, J+1); last column is the target
-    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,
-            init_weights=init_weights,
-            min_decrease=float(min_decrease),
-            max_iter=int(max_iter),
-            return_convergence=True,
-        )
+    w, converged = _sc_weight_nnls(packed, return_convergence=True)
     w = np.asarray(w, dtype=float)
     return w, y1[:n_out] - Y0[:n_out] @ w, bool(converged)
 
@@ -246,9 +234,6 @@ def _single_null_pvalue(
     perms: np.ndarray,
     q: Any,
     *,
-    max_iter: int,
-    min_decrease: float,
-    init_weights: Optional[np.ndarray] = None,
     alternative: str = "two-sided",
     x1_rows: Optional[np.ndarray] = None,
     X0_rows: Optional[np.ndarray] = None,
@@ -265,9 +250,6 @@ def _single_null_pvalue(
     w, resid, conv = _cwz_proxy_fit(
         y1n,
         Y0,
-        max_iter=max_iter,
-        min_decrease=min_decrease,
-        init_weights=init_weights,
         x1_rows=x1_rows,
         X0_rows=X0_rows,
     )
@@ -346,8 +328,6 @@ def _invert_single_post(
     alpha: float,
     perms: np.ndarray,
     *,
-    max_iter: int,
-    min_decrease: float,
     grid: Optional[np.ndarray] = None,
     n_grid: int = 100,
     alternative: str = "two-sided",
@@ -359,9 +339,8 @@ def _invert_single_post(
     The series ``y1_obs`` ``(m,)`` / ``Y0`` ``(m, J)`` has exactly one post slot at
     ``post_idx`` (a per-period CI sub-series ``(pre ∪ {t})`` or a block-collapsed
     average series). For each candidate effect ``θ`` the post outcome is imputed
-    (``y1[post_idx] − θ``), the proxy refit (warm-started across the grid — FW is
-    convex so warm-starting only affects speed, not the optimum), and the
-    permutation p-value recomputed. With one post slot ``S_q`` reduces to
+    (``y1[post_idx] − θ``), the proxy refit exactly (one NNLS per grid point; no
+    warm start is needed or used), and the permutation p-value recomputed. With one post slot ``S_q`` reduces to
     ``|û_post|`` for every ``q``, so ``q`` is inert here (fixed to 1); under
     ``alternative="greater"``/``"less"`` the single-slot statistic is the SIGNED
     ``û_post`` (resp. ``−û_post``), and the inversion yields a ONE-SIDED
@@ -391,8 +370,6 @@ def _invert_single_post(
     w_pre, resid_pre, _ = _cwz_proxy_fit(
         y1_obs[pre_sel],
         Y0[pre_sel],
-        max_iter=max_iter,
-        min_decrease=min_decrease,
         x1_rows=x1_rows,
         X0_rows=X0_rows,
     )
@@ -420,20 +397,15 @@ def _invert_single_post(
 
     pvals = np.empty(grid.shape[0], dtype=float)
     converged = np.empty(grid.shape[0], dtype=bool)
-    w_prev = w_pre
     for i, theta in enumerate(grid):
         y1n = y1_obs.copy()
         y1n[post_idx] = y1n[post_idx] - float(theta)
         w, resid, conv = _cwz_proxy_fit(
             y1n,
             Y0,
-            max_iter=max_iter,
-            min_decrease=min_decrease,
-            init_weights=w_prev,
             x1_rows=x1_rows,
             X0_rows=X0_rows,
         )
-        w_prev = w
         p, _, _ = _cwz_pvalue(resid, post_mask, perms, 1, alternative=alternative)
         pvals[i] = p
         converged[i] = conv

diff_diff/guides/llms-full.txt

@@ -645,8 +645,8 @@ SyntheticControl(
     custom_v=None,                          # diagonal V (len = #predictors); required iff v_method="custom"; forbidden otherwise
     optimizer_options: dict | None = None,  # merged into scipy.optimize.minimize (outer V search; nested/cv)
     n_starts: int = 4,                      # multistart count for the outer V search (nested/cv)
-    inner_max_iter: int = 10000,            # Frank-Wolfe inner-solve cap
-    inner_min_decrease: float = 1e-5,       # inner-solve convergence scale (scale-aware)
+    inner_max_iter: int | None = None,      # DEPRECATED no-op (warns; removed in 5.0): the inner solve is exact (NNLS)
+    inner_min_decrease: float | None = None,  # DEPRECATED no-op (warns; removed in 5.0): the inner solve is exact (NNLS)
     standardize: str = "std",               # "std" (per-row SD, ddof=1) or "none" (deviation from R)
     alpha: float = 0.05,
     seed: int | None = None,                # seeds the multistart Dirichlet draws