Worked examples

Each example is a standalone script using only skgrad’s runtime dependencies. Install the checkout, then run python scripts/check_examples.py to execute all examples and their numerical assertions. These checks verify derivative semantics, not predictive quality of the small fitted demonstration models.

Binary decision scores

The quick-start example shows a constant affine gradient. Binary classification returns the score for the positive class, not a probability derivative.

import numpy as np
from sklearn.linear_model import LogisticRegression

import skgrad


X = np.array([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]])
y = np.array([0, 1, 0, 1])
model = LogisticRegression().fit(X, y)

scores, jacobian = skgrad.value_and_jacobian(model, [[0.25, 0.75]])
print(scores)
print(jacobian)

Multiclass score gradients

A target is an output position in classes_. This multinomial logistic example checks the score outputs and their softmax relationship to sklearn probabilities. The gradients themselves remain on the score scale.

"""Differentiate each class score and verify softmax probabilities."""
import numpy as np
from sklearn.linear_model import LogisticRegression
from scipy.special import softmax
import skgrad

rng = np.random.default_rng(12)
X = rng.normal(size=(150, 3))
y = np.argmax(np.column_stack((X[:, 0], X[:, 1], -X[:, 0])), axis=1)
model = LogisticRegression(max_iter=500).fit(X, y)
values, jacobian = skgrad.value_and_jacobian(model, X[:4])
np.testing.assert_allclose(values, model.decision_function(X[:4]))
np.testing.assert_allclose(softmax(values, axis=1), model.predict_proba(X[:4]))
for target, label in enumerate(model.classes_):
    gradient = skgrad.input_gradient(model, X[:4], target=target)
    np.testing.assert_allclose(gradient, jacobian[:, target, :])
    print(f"Class {label}: first input gradient {gradient[0]}")

Multiple neural-network outputs

Use target=1 for the second regression output. Selected MLP gradients avoid forming every output’s Jacobian; here the complete result is computed only to verify equivalence. Expected shapes are (5, 2), (5, 2, 3), and (5, 3).

"""Request one MLP output without constructing all gradients."""
import numpy as np
from sklearn.neural_network import MLPRegressor
import skgrad

rng = np.random.default_rng(7)
X = rng.normal(size=(120, 3))
y = np.column_stack((np.sin(X[:, 0]), X[:, 1] * X[:, 2]))
model = MLPRegressor(hidden_layer_sizes=(12,), activation="tanh",
                     solver="lbfgs", max_iter=2000, random_state=7).fit(X, y)
values, jacobian = skgrad.value_and_jacobian(model, X[:5])
selected = skgrad.input_gradient(model, X[:5], target=1)
np.testing.assert_allclose(values, model.predict(X[:5]))
np.testing.assert_allclose(selected, jacobian[:, 1, :], atol=1e-12)
print("Values, full Jacobian, selected gradient:", values.shape,
      jacobian.shape, selected.shape)

Checking a kernel gradient

Central differences evaluate sklearn’s own prediction function, providing an independent numerical comparison. The check uses smooth RBF regression in float64, with a fixed step and explicit tolerances.

"""Check an RBF SVR gradient against sklearn's own predictions."""
import numpy as np
from sklearn.svm import SVR
import skgrad

rng = np.random.default_rng(8)
X = rng.normal(size=(80, 3))
y = np.sin(X[:, 0]) + X[:, 1] * X[:, 2]
model = SVR(kernel="rbf", gamma=0.4).fit(X, y)
evaluation = X[:4]
analytic = skgrad.input_gradient(model, evaluation)
step = 1e-5
numerical = np.empty_like(evaluation)
for feature in range(evaluation.shape[1]):
    plus, minus = evaluation.copy(), evaluation.copy()
    plus[:, feature] += step
    minus[:, feature] -= step
    numerical[:, feature] = (model.predict(plus) - model.predict(minus)) / (2 * step)
np.testing.assert_allclose(analytic, numerical, atol=1e-8, rtol=1e-6)
print("Maximum absolute error:", np.max(np.abs(analytic - numerical)))

Scaling inputs explicitly

For z_j = (x_j - mean_j) / scale_j, the original-coordinate derivative is df/dx_j = (df/dz_j) / scale_j. The pipeline now performs this chain rule automatically; the example checks its result against manual scaling and independent finite differences. It assumes StandardScaler’s default with_std=True and continuous inputs.

"""Compose a scaler explicitly and return derivatives in original units."""
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import Ridge
from sklearn.pipeline import make_pipeline
import skgrad

rng = np.random.default_rng(9)
X = rng.normal(size=(100, 2)) * np.array([2.0, 20.0])
y = 3 * X[:, 0] - 0.5 * X[:, 1]
scaler = StandardScaler().fit(X)
model = Ridge().fit(scaler.transform(X), y)
evaluation = X[:3]
scaled_gradient = skgrad.input_gradient(model, scaler.transform(evaluation))
original_gradient = scaled_gradient / scaler.scale_[None, :]
pipeline = make_pipeline(scaler, model)
assert skgrad.supports(pipeline)
np.testing.assert_allclose(
    skgrad.input_gradient(pipeline, evaluation), original_gradient, atol=1e-12
)
step = 1e-4
for feature in range(2):
    plus, minus = evaluation.copy(), evaluation.copy()
    plus[:, feature] += step
    minus[:, feature] -= step
    numerical = (model.predict(scaler.transform(plus))
                 - model.predict(scaler.transform(minus))) / (2 * step)
    np.testing.assert_allclose(original_gradient[:, feature], numerical, atol=1e-8)
print("Gradient in original input units:", original_gradient[0])

Integrating a polynomial gradient

A cubic polynomial has quadratic derivatives along a straight path. Two Gauss–Legendre points integrate those derivatives exactly up to floating-point error. Post-expansion scaling is already included by skgrad. Multiplication by the input displacement gives feature contributions whose sum equals the fitted prediction difference. The zero baseline here is illustrative, not a universal choice of meaningful reference.

"""Integrate an analytic polynomial gradient along one straight path."""
import numpy as np
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.linear_model import Ridge
import skgrad

rng = np.random.default_rng(10)
X = rng.normal(size=(150, 2))
y = X[:, 0] ** 3 + X[:, 0] * X[:, 1] - 2 * X[:, 1]
model = make_pipeline(PolynomialFeatures(3), StandardScaler(), Ridge(alpha=0.01))
model.fit(X, y)
baseline = np.zeros(2)
point = np.array([0.8, -0.4])
order = skgrad.gradient_properties(model).exact_quadrature_steps
nodes, weights = np.polynomial.legendre.leggauss(order)
path = baseline + ((nodes + 1) / 2)[:, None] * (point - baseline)
gradients = skgrad.input_gradient(model, path)
attributions = (point - baseline) * ((weights / 2) @ gradients)
difference = model.predict(point[None, :])[0] - model.predict(baseline[None, :])[0]
np.testing.assert_allclose(attributions.sum(), difference, atol=1e-12)
print("Quadrature points:", order)
print("Feature contributions:", attributions)
print("Completeness residual:", attributions.sum() - difference)

Scaling and PCA before an MLP

This nested pipeline reduces four features to three PCA components. Its returned gradient still has four columns in original input order. The script independently checks those derivatives against perturbations of the complete sklearn pipeline.

"""Differentiate a complete scaled/PCA/MLP pipeline in original coordinates."""
import numpy as np
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.neural_network import MLPRegressor
import skgrad

rng = np.random.default_rng(17)
X = rng.normal(size=(100, 4))
y = np.sin(X[:, 0]) + X[:, 1] - X[:, 2]
preprocessing = make_pipeline(StandardScaler(), PCA(3, whiten=True))
model = make_pipeline(preprocessing, MLPRegressor(
    hidden_layer_sizes=(8,), activation="tanh", solver="lbfgs",
    max_iter=2000, tol=1e-3, alpha=0.1, random_state=17,
)).fit(X, y)
evaluation = X[:3]
gradient = skgrad.input_gradient(model, evaluation)
assert gradient.shape == (3, 4)  # original inputs, not three PCA components
step = 1e-5
for feature in range(4):
    plus, minus = evaluation.copy(), evaluation.copy()
    plus[:, feature] += step
    minus[:, feature] -= step
    numerical = (model.predict(plus) - model.predict(minus)) / (2 * step)
    np.testing.assert_allclose(gradient[:, feature], numerical, atol=1e-7)
print("Gradient in original feature coordinates:", gradient)

Choosing original or standardized gradients

Both views accept original inputs, but return derivatives in the explicitly selected coordinates. Their output values agree; their gradient units differ.

"""Compare original and standardized gradients with one explicit boundary."""
import numpy as np
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import Ridge
import skgrad

rng = np.random.default_rng(11)
X = rng.normal(size=(60, 3)) * [1, 4, 9]
model = Pipeline([("scale", StandardScaler()), ("regressor", Ridge())]).fit(
    X, X[:, 0] - X[:, 1]
)
original = skgrad.pipeline_view(model)
standardized = skgrad.pipeline_view(model, after="scale")
# Both view methods accept the same ORIGINAL inputs.
raw_gradient = original.input_gradient(X[:4])
z_gradient = standardized.input_gradient(X[:4])
np.testing.assert_allclose(z_gradient / model[0].scale_, raw_gradient)
np.testing.assert_allclose(original.model_output(X[:4]), standardized.model_output(X[:4]))
print("Original gradient:", raw_gradient[0])
print("Standardized gradient:", z_gradient[0])
print("Selected names:", standardized.get_feature_names_out())