The v6.0 pipeline ranks 27 features by SHAP importance, which tells you what matters. This post is about getting an explicit, written-on-paper equation instead — what matters and how.
Post 8 closed with a stacking ensemble, SHAP attributions, and a bootstrap-calibrated confidence interval on every prediction. That's a complete black-box pipeline. But SHAP importance answers “which of the 27 features mattered for this prediction,” not “what is the relationship.” You can't write a SHAP plot on a whiteboard, compare its coefficients to the GKA prefactor, or check whether it respects the Mott criterion outside the training range. Post 1's roadmap called this step “looking for a formula,” and this is the first real attempt at it.
Why SHAP isn't enough on its own
SHAP values are local and model-specific — they explain one prediction from one trained random forest, and a different forest, or a different random seed, can redistribute importance among correlated features without changing the predictions much. None of that is a flaw in the v6.0 pipeline; it's just a different question than the one this post is asking. A formula is a single, fixed functional form with fixed coefficients that you can inspect, falsify, and compare against the GKA/Mott/ZSA limits the whole series has been building toward. The cost is that finding one is harder and easier to fool yourself with than ranking feature importances.
Picking a method
Three reasonable options, in order of how much they cost to run:
| Method | What it does | Trade-off |
|---|---|---|
| Sparse regression (SISSO-style) | Build a library of physics-motivated candidate terms (products, ratios, powers of existing features), then use L1-regularized regression to pick a sparse subset. | Fully transparent, runs anywhere sklearn does, no new dependency. You have to propose the candidate terms yourself. |
| gplearn | Genetic-programming symbolic regression in pure Python/sklearn style — evolves expression trees from a function set. | Finds nonlinear forms you wouldn't think to propose, but the result needs more scrutiny — expressions can get baroque before parsimony pruning catches up. |
| PySR | State-of-the-art genetic-programming symbolic regression, generally the best accuracy/simplicity frontier. | Needs a Julia backend underneath — extra friction on Colab/Kaggle, where this whole pipeline already auto-detects platform. |
Given that the pipeline already runs across local Python, Colab, and Kaggle, I'm starting with sparse regression: zero new dependencies, and it lets the candidate library be biased toward the physics this series already cares about — including the ΔCT×angle interaction terms post 10 flagged as the coupling test. gplearn is the fallback if the sparse library turns up nothing convincing, run as a sanity check on forms outside the candidate library, not as the primary method.
Application: building the candidate library and fitting it
This isn't a toy widget like posts 5, 6, 9, and 10 — it's the actual code, meant to run against the real feature table the v6.0 pipeline already exports. Two steps: build the candidate library, then let an information-criterion-regularized regression pick the sparse subset.
# from the 27 features the v6.0 pipeline already exports (Cell 9)
import pandas as pd
from itertools import combinations
from sklearn.linear_model import LassoLarsIC
from sklearn.preprocessing import StandardScaler
df = pd.read_csv("mxcy_features_v6.csv")
target = "band_gap_eV" # or "magnetization_uB"
exclude = {target, "compound", "ordering_label"}
base_features = [c for c in df.columns if c not in exclude]
# candidate library: every base feature, plus pairwise products and
# ratios — this is where post 10's coupling terms (Δ_CT × angle, etc.)
# get the chance to show up if they actually matter
candidates = {f: df[f] for f in base_features}
for f1, f2 in combinations(base_features, 2):
candidates[ff"{f1}*{f2}"] = df[f1] * df[f2]
if (df[f2] != 0).all():
candidates[ff"{f1}/{f2}"] = df[f1] / df[f2]
X = pd.DataFrame(candidates)
y = df[target]
X_scaled = StandardScaler().fit_transform(X) # coefficients become comparable
# this is what actually makes the result a FORMULA, not a ranking —
# most coefficients get driven exactly to zero
model = LassoLarsIC(criterion="bic")
model.fit(X_scaled, y)
terms = [(name, coef) for name, coef in zip(X.columns, model.coef_)
if abs(coef) > 1e-6]
terms.sort(key=lambda t: -abs(t[1]))
print(ff"{len(terms)} surviving terms, intercept={model.intercept_:.4f}")
for name, coef in terms:
print(ff" {coef:+.4f} * {name}")
LassoLarsIC is doing the actual symbolic-regression work here: the L1 penalty plus BIC model-order selection drives most of the hundreds of candidate terms exactly to zero, leaving a short, inspectable sum. That's the formula — not the full candidate library, just whatever survives.
Reading a result without fooling yourself
A library this size relative to a few hundred compounds is exactly the setup where it's easy to find a formula that fits the training set and means nothing. Three checks before trusting any surviving term: hold out a test split (or better, cross-validate the BIC selection itself, not just the final coefficients) and confirm the same terms keep surviving; check the sign of each coefficient against the physics it's supposed to represent — if the ΔCT×angle term comes out the wrong sign relative to what post 10 predicted, that's a flag for a confound or leakage, not a discovery; and watch for ratio terms with a denominator that's near zero anywhere in the dataset, since those can produce huge, unstable coefficients that are fitting artifacts of one or two compounds rather than real structure.
What happens next
Whatever formula survives those checks isn't the end point — it's the input to the last post in this arc: validating it explicitly against the GKA, Mott-Hubbard, and ZSA limits this entire series has been building descriptors around, and compiling whatever holds up into a single interpretable rule bank.
My name is Dr Abderrahmane Reggad. I'm a 54 year old and I am a researcher in materials' sciences.
0 Comments