Training-Serving Skew and Feature Stores

Training-serving skew occurs when the inputs a model receives in service differ from the inputs used to train or evaluate it because of their construction, timing, or representation. A historical aggregate that includes future outcomes, a stale user profile, and a unit conversion applied only at serving time are different causes. A change in the population can also change performance even when the feature code is identical. The examples below separate these possibilities using simulations, not a deployed service.

Run the body examples in order. They share the simulated events and history arrays. Model fitting uses the earlier 70% of events; evaluation uses the later 30%. Early stopping is disabled so the demonstration does not introduce an additional random internal validation split.

A feature computed as of when?

A historical conversion rate needs a cutoff and a definition of which outcomes were available by then. The simulation orders rows by prediction time and assumes each label becomes available immediately after its prediction. Thus later rows may use earlier test-period outcomes, but never their own or future labels. For delayed conversions, an availability-time cutoff is needed as well. The smoothed rate is \((s+20\times0.5)/(c+20)\), where s is the number of prior successes and c the number of prior requests. A new user starts at 0.5; after one success the value is 11/21. The prior 0.5 is fixed in advance, not estimated from the full dataset.

import numpy as np
from collections import defaultdict
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.metrics import roc_auc_score

g = np.random.default_rng(0)
n_users, n_ev = 4000, 40000
uid = g.integers(0, n_users, n_ev)
user_q = g.normal(size=n_users)                      # latent per-user quality
x = g.normal(size=(n_ev, 3))
y = (g.random(n_ev) < 1 / (1 + np.exp(-(0.9 * user_q[uid] + 0.7 * x[:, 0]
                                        - 0.3 * x[:, 1])))).astype(int)

prior_m, prior_n = 0.5, 20.0
past = np.zeros(n_ev); tot_s = np.zeros(n_users); tot_c = np.zeros(n_users)
seq = defaultdict(list)
for i in range(n_ev):                                # outcome is available immediately after this row
    u = uid[i]
    past[i] = (tot_s[u] + prior_m * prior_n) / (tot_c[u] + prior_n)
    seq[u].append(i); tot_s[u] += y[i]; tot_c[u] += 1

gs = np.zeros(n_users); gc = np.zeros(n_users)       # WHOLE-HISTORY: what a GROUP BY gives
np.add.at(gs, uid, y); np.add.at(gc, uid, 1)
full = (gs[uid] + prior_m * prior_n) / (gc[uid] + prior_n)
split = int(0.7 * n_ev)

def fit_eval(train_feat, serve_feat, label):
    m = HistGradientBoostingClassifier(random_state=0, early_stopping=False).fit(
        np.c_[x[:split], train_feat[:split]], y[:split])
    off = roc_auc_score(y[split:], m.predict_proba(np.c_[x[split:], train_feat[split:]])[:, 1])
    on = roc_auc_score(y[split:], m.predict_proba(np.c_[x[split:], serve_feat[split:]])[:, 1])
    print(f"{label:34s} validation-input AUC {off:.4f}   serving-input replay AUC {on:.4f}   gap {off - on:+.4f}")

fit_eval(full, past, "trained on whole-history feature")
fit_eval(past, past, "trained on past-only feature")
# trained on whole-history feature   validation-input AUC 0.8067   serving-input replay AUC 0.7161   gap +0.0906
# trained on past-only feature       validation-input AUC 0.7195   serving-input replay AUC 0.7195   gap +0.0000

The whole-history feature uses every event, including each row’s own label and later labels. The first model scores 0.8067 with that feature and 0.7161 with the available-history input, a gap of 0.0906. Holding out model-fitting rows does not undo this feature leakage. The serving-input column recomputes scores on the same held-out labels using only the simulated available history. This replay exposes a discrepancy that the leaky validation inputs conceal; a properly time-restricted feature pipeline or an audit of feature construction can detect it before deployment.

The second model is trained and evaluated with the same past-only construction. Its two reported scores are equal because the arrays supplied to the model are identical. That equality is an implementation check in this replay, not a guarantee of future production performance. Comparing the two models on available-history inputs also measures the effect of training on the unavailable feature in this particular run; it does not establish a universal degradation size.

Logging the feature vector actually supplied to a model provides a useful replay source, but logging alone does not guarantee correctness. Missing events, delayed labels, identity mismatches, version changes, and filtering can still produce a biased training set. Shared transformations, time-correct historical reconstruction, and comparisons against served values are complementary controls. A logged feature can itself already contain leakage or the wrong units.

What a feature store must keep consistent

A feature definition needs an entity key, value type and units, aggregation window, missing-value policy, and transformation version. An offline store retains history for training and replay; an online store serves values with low latency. A registry records definitions and metadata. Batch materialization copies computed feature values into the serving store; streaming jobs can update values as events arrive. Shared definitions help consistency, but the system still needs checks on data arrival, computation, and retrieval. A feature store is not a guarantee against leakage.

A historical join must use the right entity and prediction cutoff. Event time says when the source event happened; availability time says when the pipeline could actually use it. A conversion occurring at hour 4 but arriving at hour 8 must not appear in a replay of a prediction at hour 6. Later corrections and backfills also need versioned availability information. Selecting only by event time can reconstruct a history more complete than the service actually had.

The records below are versioned feature snapshots. as_of is the snapshot’s source cutoff; available_at is when that snapshot became usable. Times are hours from a common origin. A snapshot is eligible only if both times are at or before the prediction time; ties between availability and prediction require an explicit ordering rule in a real system. Among eligible snapshots of the requested version, this toy lookup selects the latest source cutoff, then the latest arrival for that cutoff.

snapshots = [
    dict(entity="A", as_of=1, available_at=2, version=1, value=2),
    dict(entity="A", as_of=4, available_at=8, version=1, value=5),
    dict(entity="A", as_of=10, available_at=10, version=1, value=7),
    dict(entity="A", as_of=4, available_at=5, version=2, value=50),
]

def lookup(records, entity, prediction_time, version):
    candidates = [row for row in records
                  if row["entity"] == entity and row["version"] == version
                  and row["as_of"] <= prediction_time
                  and row["available_at"] <= prediction_time]
    if not candidates:
        return None
    return max(candidates, key=lambda row: (row["as_of"], row["available_at"]))["value"]

for time in (6, 9):
    print(f"A at hour {time}, version 1: {lookup(snapshots, 'A', time, 1)}")
print("A at hour 6, version 2:", lookup(snapshots, "A", 6, 2))
print("unknown entity:", lookup(snapshots, "B", 9, 1))
# A at hour 6, version 1: 2
# A at hour 9, version 1: 5
# A at hour 6, version 2: 50
# unknown entity: None

The hour-6 lookup returns 2 even though an event-time-only lookup could return 5. At hour 9, version 1 returns 5. Version 2 is a separate definition and must not be silently substituted. None means unavailable, not a measured zero. Production retrieval needs a declared age limit and fallback policy, plus monitoring for default use. A time-to-live limit (TTL) used in historical retrieval is not automatically a guarantee of online freshness; check the store’s semantics.

This scan illustrates the selection rules; it is not a feature-store implementation. A historical store might reconstruct snapshots from immutable event and arrival logs, while an online store may retain only the latest value. Replay must model the actual update, caching, expiry, and fallback behavior if the goal is to reproduce served inputs. Log request identifiers, feature and model versions, source cutoffs, retrieval times, and fallback reasons where needed to make that comparison possible.

How fresh does a feature have to be?

m = HistGradientBoostingClassifier(random_state=0, early_stopping=False).fit(
    np.c_[x[:split], past[:split]], y[:split])
print(f"{'staleness (user events)':>24} {'replay AUC':>11} {'prior-only fraction':>20}")
for k in (0, 1, 3, 10, 30, 100):
    stale = np.empty(n_ev)
    initial = np.empty(n_ev, dtype=bool)
    for u, idxs in seq.items():                       # the value from k of this user's events ago
        arr = past[idxs]
        historical_index = np.maximum(np.arange(len(idxs)) - k, 0)
        stale[idxs] = arr[historical_index]
        initial[idxs] = historical_index == 0
    auc = roc_auc_score(y[split:], m.predict_proba(np.c_[x[split:], stale[split:]])[:, 1])
    print(f"{k:24d} {auc:11.4f} {initial[split:].mean():20.4f}")
#  staleness (user events)  replay AUC  prior-only fraction
#                        0      0.7195               0.0004
#                        1      0.7159               0.0022
#                        3      0.7061               0.0348
#                       10      0.6753               0.7590
#                       30      0.6681               1.0000
#                      100      0.6681               1.0000

The table replays one fixed model using values from k earlier events for each user. It measures event-count lag, not seconds or hours. The prior-only fraction shows how much of the held-out population has been reduced to its initial value. At large enough k, all values become the same prior; the resulting score reflects how this fitted model behaves when that feature is constant.

The fresh-versus-prior gap is specific to this fitted model, dataset, and perturbation. It is not the feature’s total value, the score of a model retrained without it, or an upper bound on every future staleness cost. To choose a refresh interval, measure actual feature ages and replay realistic delays, accounting for user activity and delayed events. Compare quality and downstream costs with the latency and operating cost of the proposed pipeline.

The bug that does not raise an exception

from sklearn.linear_model import LogisticRegression

m2 = LogisticRegression(max_iter=2000).fit(np.c_[x[:split], past[:split]], y[:split])
base_scores = m2.decision_function(np.c_[x[split:], past[split:]])
base = roc_auc_score(y[split:], base_scores)
for mult in (1.0, 0.001, 1000.0, -1.0):
    xs = x[split:].copy()
    xs[:, 0] *= mult
    scores = m2.decision_function(np.c_[xs, past[split:]])
    auc = roc_auc_score(y[split:], scores)
    loss = np.mean(np.logaddexp(0, scores) - y[split:] * scores)
    print(f"x0 multiplier {mult:8g}  score AUC {auc:.4f}  AUC drop {base-auc:+.4f}  log loss {loss:.4f}")
# x0 multiplier        1  score AUC 0.7285  AUC drop +0.0000  log loss 0.6081
# x0 multiplier    0.001  score AUC 0.6589  AUC drop +0.0696  log loss 0.6548
# x0 multiplier     1000  score AUC 0.6611  AUC drop +0.0674  log loss 162.8531
# x0 multiplier       -1  score AUC 0.5120  AUC drop +0.2165  log loss 0.7885

These inputs have no physical unit; multiplying the first feature by 0.001 or 1000 mimics an inconsistent unit conversion. The code reports AUC from the linear decision score and a numerically stable logistic loss. Using the decision score for AUC avoids artificial ties from probabilities rounding to 0 or 1. At a multiplier of 1000, score AUC is 0.6611 and log loss is 162.8531. AUC measures ranking, while log loss also reflects probability quality; extreme score scaling can damage the latter substantially. Compatible array shapes and numeric types allow these examples to execute, although range or distribution checks may detect the altered values.

Exercises

1. Feature parity and model impact. Introduce six serving-input changes and compare the fraction of differing rows with each change’s measured effect on AUC.

Nearly all values can differ while the measured AUC changes very little. Difference frequency does not measure severity.

import numpy as np
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.metrics import roc_auc_score

g = np.random.default_rng(0)
n, d = 40000, 8
X = g.normal(size=(n, d))
w = np.array([1.1, -0.9, 0.6, 0.4, 0.25, 0.15, 0.05, 0.0])
y = (g.random(n) < 1 / (1 + np.exp(-(X @ w)))).astype(int)
split = int(0.7 * n)
m = HistGradientBoostingClassifier(random_state=0, early_stopping=False).fit(X[:split], y[:split])
Xte = X[split:]
base = roc_auc_score(y[split:], m.predict_proba(Xte)[:, 1])

bugs = [("2% scale drift on x2", 2, lambda v: v * 1.02),
        ("2% scale drift on x0", 0, lambda v: v * 1.02),
        ("rounded to 1 dp at serving (x5)", 5, lambda v: np.round(v, 1)),
        ("offset +0.5 on x0", 0, lambda v: v + 0.5),
        ("sign flipped on x0", 0, lambda v: -v),
        ("column dropped -> NaN (x7)", 7, lambda v: np.full_like(v, np.nan))]
print(f"{'bug':34s} {'rows differ':>12} {'AUC':>8} {'damage':>8}")
for name, j, f in bugs:
    s = Xte.copy(); s[:, j] = f(Xte[:, j])
    diff = ~np.isclose(Xte[:, j], s[:, j], rtol=1e-5, atol=1e-8, equal_nan=True)
    auc = roc_auc_score(y[split:], m.predict_proba(s)[:, 1])
    print(f"{name:34s} {diff.mean():12.4f} {auc:8.4f} {base - auc:9.6f}")
# bug                                 rows differ      AUC   damage
# 2% scale drift on x2                     1.0000   0.8166  0.000104
# 2% scale drift on x0                     1.0000   0.8169 -0.000142
# rounded to 1 dp at serving (x5)          0.9998   0.8167  0.000057
# offset +0.5 on x0                        1.0000   0.8160  0.000766
# sign flipped on x0                       1.0000   0.5178  0.298886
# column dropped -> NaN (x7)               1.0000   0.8172 -0.000505

2. Sensitivity to constant replacement. Replace each feature with its training-set mean while keeping the model fixed. Rank the resulting AUC drops. What does this perturbation say about monitoring priorities?

This ranking describes one failure scenario for this fitted model; it is not a complete ranking of feature importance or operational risk.

import numpy as np
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.metrics import roc_auc_score

g = np.random.default_rng(0)
n, d = 40000, 8
X = g.normal(size=(n, d))
w = np.array([1.1, -0.9, 0.6, 0.4, 0.25, 0.15, 0.05, 0.0])
y = (g.random(n) < 1 / (1 + np.exp(-(X @ w)))).astype(int)
split = int(0.7 * n)
m = HistGradientBoostingClassifier(random_state=0, early_stopping=False).fit(X[:split], y[:split])
Xte = X[split:]
base = roc_auc_score(y[split:], m.predict_proba(Xte)[:, 1])

print(f"{'column':>7} {'|w|':>7} {'AUC if replaced by its mean':>29} {'damage':>9}")
rows = []
for j in range(d):
    one = Xte.copy(); one[:, j] = X[:split, j].mean()
    auc = roc_auc_score(y[split:], m.predict_proba(one)[:, 1])
    rows.append((base - auc, j, auc))
for dmg, j, auc in sorted(rows, reverse=True):
    print(f"{j:7d} {abs(w[j]):7.2f} {auc:29.4f} {dmg:9.4f}")
#  column     |w|   AUC if replaced by its mean    damage
#       0    1.10                        0.7250    0.0917
#       1    0.90                        0.7608    0.0559
#       2    0.60                        0.7911    0.0256
#       3    0.40                        0.8056    0.0111
#       4    0.25                        0.8115    0.0052
#       5    0.15                        0.8168   -0.0000
#       7    0.00                        0.8173   -0.0005
#       6    0.05                        0.8173   -0.0006

3. Why the gap does not add up. Apply three separate causes of an offline-online gap one at a time and then together.

Use the same noise realization in the individual and combined scenarios, then measure the offset’s contribution with and without the other changes.

import numpy as np
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.metrics import roc_auc_score

g = np.random.default_rng(0)
n, d = 60000, 6
X = g.normal(size=(n, d))
w = np.array([1.1, -0.9, 0.6, 0.4, 0.25, 0.0])
y = (g.random(n) < 1 / (1 + np.exp(-(X @ w)))).astype(int)
split = int(0.6 * n)
m = HistGradientBoostingClassifier(random_state=0, early_stopping=False).fit(X[:split], y[:split])
Xte, yte = X[split:], y[split:]
offline = roc_auc_score(yte, m.predict_proba(Xte)[:, 1])
print(f"offline AUC (clean holdout)  {offline:.4f}")

skew = Xte.copy(); skew[:, 0] += 0.5                        # a serving offset bug
sel = np.zeros(len(Xte), bool)
sel[np.argsort(Xte[:, 1])[: len(Xte) // 2]] = True          # population shift
noise = g.normal(0, 0.35, Xte.shape)
noisier = Xte + noise                # noisier live inputs

a_skew = roc_auc_score(yte, m.predict_proba(skew)[:, 1])
a_shift = roc_auc_score(yte[sel], m.predict_proba(Xte[sel])[:, 1])
a_noise = roc_auc_score(yte, m.predict_proba(noisier)[:, 1])
both = skew + noise
a_without_offset = roc_auc_score(yte[sel], m.predict_proba(noisier[sel])[:, 1])
a_all = roc_auc_score(yte[sel], m.predict_proba(both[sel])[:, 1])

for lbl, v in (("serving offset bug only", a_skew), ("population shift only", a_shift),
               ("noisier live inputs only", a_noise), ("shift and noise, no offset", a_without_offset), ("all three together", a_all)):
    print(f"  {lbl:28s} {v:.4f}   ({v - offline:+.4f})")
print(f"  {'sum of individual effects':28s}          "
      f"({(a_skew - offline) + (a_shift - offline) + (a_noise - offline):+.4f})")
print(f"offset contribution with shift and noise {a_all - a_without_offset:+.4f}")
# offline AUC (clean holdout)  0.8199
#   serving offset bug only      0.8196   (-0.0003)
#   population shift only        0.7940   (-0.0259)
#   noisier live inputs only     0.7982   (-0.0217)
#   shift and noise, no offset   0.7698   (-0.0501)
#   all three together           0.7691   (-0.0508)
#   sum of individual effects             (-0.0479)
# offset contribution with shift and noise -0.0007

References

Feast components, feature definitions, and point-in-time joins describe one implementation of these ideas. Availability-aware replay additionally requires preserving the information about when data became usable; do not assume an event timestamp alone provides it.