gcm Shapley least-squares estimation broken on scikit-learn 1.9 (extreme sample weights)

[immu4989]

Description

The EXACT_FAST and SUBSET_SAMPLING Shapley approximation methods in dowhy.gcm.shapley return incorrect values on scikit-learn 1.9.0. Three tests fail:

• test_given_many_features_when_estimate_shapley_values_exact_fast_then_returns_correct_linear_shapley_values
• test_given_many_features_when_estimate_shapley_values_with_subset_sampling_then_returns_correct_linear_shapley_values
• test_given_specific_random_seed_when_estimate_shapley_values_with_subset_sampling_then_returns_deterministic_result

Caught by the latest-deps CI work in #1665 / #1672 (fails there on every run with scikit-learn 1.9.0; passes on 1.8.0).

Root cause

_approximate_shapley_values_via_least_squares_regression solves a weighted least squares problem via LinearRegression().fit(..., sample_weight=weights) where the full/empty subsets get weight 10**20 (to approximate equality constraints) and all other subsets weight ~1. scikit-learn 1.9 changed LinearRegression's solve path for sample weights, and with weights spanning 20 orders of magnitude the solution now loses all precision.

Minimal repro independent of DoWhy:

import numpy as np
from sklearn.linear_model import LinearRegression
rng = np.random.RandomState(0)
X = (rng.uniform(size=(500, 12)) < 0.5).astype(float)
X[0] = 0.0; X[1] = 1.0
beta = rng.normal(size=12)
y = X @ beta + 3.0
w = np.ones(500); w[0] = w[1] = 1e20
coef = LinearRegression().fit(X, y, sample_weight=w).coef_
print(np.max(np.abs(coef - beta)))
# scikit-learn 1.8.0: ~1e-14
# scikit-learn 1.9.0: ~2.08

Fix

Solve the weighted least squares explicitly (weighted centering + np.linalg.lstsq on the square-root-weighted system), which reproduces the numerically stable behavior and makes the result independent of solver changes in third-party libraries. PR: #1699.

Note

test_anomaly_attribution.py::...feature_relevance... also fails in the latest-deps CI runs, but that turns out to be a pre-existing order/RNG-sensitive test (it fails on main with scikit-learn 1.8 too depending on the random state; the @flaky retries can't save it in CI because the suite ordering makes the failure deterministic). That's a separate test-robustness issue.

[github-actions]

🤖 This is an automated response from Repo Assist.

Great find, @immu4989 — and thank you for the well-structured root-cause analysis and minimal repro!

The explanation is spot-on: scikit-learn 1.9 changed the internal solve path for LinearRegression with sample weights, and with weights spanning 20 orders of magnitude (1 vs 1020) the system was already numerically ill-conditioned in a way that 1.8's solver happened to tolerate. The explicit WLS implementation in PR #1699 is the right fix: weighted centering + np.linalg.lstsq on the square-root-weighted system is both correct and independent of any third-party solver.

One thing I'd double-check in the PR: the set_function_results.ndim > 1 branch (multi-output case, where coefficients.T is returned) — the original LinearRegression.coef_ for multi-output is already shape (n_outputs, n_features), so returning coefficients.T matches that convention. Looks right, but worth a quick unit test if it isn't already covered.

On the separate test_anomaly_attribution::feature_relevance flakiness you noted: I've opened a companion PR that adds np.random.seed(0) to that specific test. Since all randomness in the test path flows through numpy's global RNG (data generation, BayesianGaussianMixture EM, permute_features), seeding once at the top of the test makes it independent of CI test ordering, which is the root cause of the @flaky-proof failures you observed.

Generated by 🌈 Repo Assist, see workflow run. Learn more.

To install this agentic workflow, run

gh aw add githubnext/agentics/workflows/repo-assist.md@11c9a2c442e519ff2b427bf58679f5a525353f76