SupraLocal Interpretable Model Agnostic Explanations (SLIME)
In this demo we explain the SLIME algorithm and show you how to:
Use
SymbolicModel’s SLIME functionalityApproximate local behavior of a black-box model
Make predictions with the surrogate symbolic model
SLIME is a very similar model interpretability technique to LIME (Local Interpretable Model Agnostic Explanations). Essentially, LIME is a technique that approximates any black box machine learning model with a linear model around a certain point. The surrogate model is fitted with Gaussian perturbed samples around the point.
SLIME extends this idea. Here, we implement the framework as presented by the SLIME paper (Fong, Motani 2025).
SLIME Algorithm: Suppose we have a black-box model \(f(x)\),
We want to fit a local surrogate model \(g^{*}(x) \in \mathcal{G}\) where \(\mathcal{S}\) is the space of closed-form analytic expressions (symbolic models) to estimate the behaviour of \(f(x)\) around point \(x_0\).
To choose our surrogate model, we create a dataset points near \(x_0\). This dataset, \(\mathcal{D}\), consists of the \(J-\) nearest neighbours around \(x_0\) and \(N_{\text{synthetic}}\) number of Gaussian sampled points around \(x_0\) (the default variance of the Gaussian is equal to half the variance of the set of \(J\) neighbours, hence we choose the size of the neighbourhood by setting \(J\)).
We evaluate the outputs of the model on this dataset, \(\mathcal{D}\), to create the targets for our symbolic regression \(\mathcal{T}\).
We choose \(g^{*}(x)\) by
where \(\pi_x(z_i) = \exp(\frac{-||x_0-x_i||^2}{\sigma^2})\) is the proximity kernel weighting and \(M\) is a weighting for the real neighbours.
SLIME has benefits over LIME as it does not restrict itself to linear models, yet does not have any reduction in interpretability of the surrogate model.
SymTorch’s SymbolicModel can fit a SLIME model for you. The user can specify:
\(J-\) nearest neighbours (
J_nn = 10)\(M\) (
real_weighting = 1.0)\(N_{\text{synthetic}}\) (
num_synthetic = 0)
As well as
Nearest neighbour metric (
nn_metric = "euclidean")PySR parameters (
pysr_params=Nonedefaults to default parameters)
The model is fitted using the .distill(SLIME = True) method and predictions can be made using .predict().
Conceptual example
Let us see an example of how SLIME can be effective at analysing black box models by using a spiral dataset.
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(7)
N = 400
theta = np.sqrt(np.random.rand(N)) * 2 * np.pi
# Class A spiral
r_a = 2 * theta + np.pi
data_a = np.column_stack([np.cos(theta) * r_a, np.sin(theta) * r_a])
x_a = data_a + np.random.randn(N, 2)
y_a = np.zeros(N, dtype=int)
# Class B spiral
r_b = -2 * theta - np.pi
data_b = np.column_stack([np.cos(theta) * r_b, np.sin(theta) * r_b])
x_b = data_b + np.random.randn(N, 2)
y_b = np.ones(N, dtype=int)
# Combine
X = np.vstack([x_a, x_b])
y = np.concatenate([y_a, y_b])
# Shuffle together
idx = np.random.permutation(len(y))
X = X[idx]
y = y[idx]
The dataset
plt.scatter(X[y==0, 0], X[y==0, 1], s=8, label="Class 0")
plt.scatter(X[y==1, 0], X[y==1, 1], s=8, label="Class 1")
plt.xlabel(r"$x_1$")
plt.ylabel(r"$x_2$")
plt.legend()
plt.title('Original dataset')
plt.show()
Training the model
Set up a simple two layer MLP.
import torch
import torch.nn as nn
import torch.nn.functional as F
class SimpleMLP(nn.Module):
def __init__(self, in_dim=2, hidden_dim=64, out_dim=2):
super().__init__()
self.layer1 = nn.Linear(in_dim, hidden_dim)
self.layer2 = nn.Linear(hidden_dim, hidden_dim)
self.output = nn.Linear(hidden_dim, out_dim)
def forward(self, x):
x = F.relu(self.layer1(x))
x = F.relu(self.layer2(x))
x = self.output(x)
return x
from torch.utils.data import TensorDataset, DataLoader
X_tensor = torch.from_numpy(X).float()
y_tensor = torch.from_numpy(y).long()
dataset = TensorDataset(X_tensor, y_tensor)
loader = DataLoader(dataset, batch_size=64, shuffle=True)
Set up training.
device = torch.device("mps")
model = SimpleMLP(in_dim=2, hidden_dim=64, out_dim=2).to(device)
criterion = nn.CrossEntropyLoss()
optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)
Train the model.
for epoch in range(50):
model.train()
running_loss = 0.0
for xb, yb in loader:
xb = xb.to(device)
yb = yb.to(device)
optimizer.zero_grad()
logits = model(xb)
loss = criterion(logits, yb)
loss.backward()
optimizer.step()
running_loss += loss.item() * xb.size(0)
epoch_loss = running_loss / len(dataset)
if (epoch + 1) % 10 == 0:
print(f"Epoch {epoch+1:02d}, loss = {epoch_loss:.4f}")
Epoch 10, loss = 0.3800
Epoch 20, loss = 0.2161
Epoch 30, loss = 0.0862
Epoch 40, loss = 0.0434
Epoch 50, loss = 0.0265
Get predictions for the data.
def model_func(inputs):
model.eval()
with torch.no_grad():
# convert numpy arrays to torch tensors
if isinstance(inputs, np.ndarray):
inputs = torch.from_numpy(inputs).float()
# ensure tensor + move to device
inputs = inputs.to(device)
logits = model(inputs)
probs = torch.softmax(logits, dim=-1)
preds = probs.argmax(dim=-1).cpu().numpy()
return preds
preds = model_func(X_tensor.to(device))
from sklearn.metrics import accuracy_score
acc = accuracy_score(y, preds)
print(f"The model accuracy on this dataset is {acc}.")
The model accuracy on this dataset is 0.9975.
Interpreting the model outputs with SLIME
Imagining that we don’t know the true form of the data, we would want to understand how the model behaves. For complex models, it is easier to see how the model behaves around a local point of interest.
np.random.seed(7)
idx = np.random.randint(len(X))
x0 = X[idx]
plt.scatter(X[y==0, 0], X[y==0, 1], s=8, label="Class 0")
plt.scatter(X[y==1, 0], X[y==1, 1], s=8, label="Class 1")
# plot the special point on top
plt.scatter(x0[0], x0[1], color="red", s=40, label=r"$x_0$")
plt.xlabel(r"$x_1$")
plt.ylabel(r"$x_2$")
plt.title(r"Original dataset with chosen local point, $x_0$")
plt.legend()
<matplotlib.legend.Legend at 0x318175dd0>
Let us just restrict ourself to the neighbourhood of points close to \(x_0\).
from sklearn.neighbors import NearestNeighbors
nbrs = NearestNeighbors(n_neighbors=50, metric='euclidean').fit(X)
_, indices = nbrs.kneighbors(x0.reshape(1, -1))
neigh_idx = indices[0]
# plot neighbours using the same class-specific colours
plt.scatter(X[neigh_idx][y[neigh_idx] == 0, 0],
X[neigh_idx][y[neigh_idx] == 0, 1],
s=20, label="Neighbours (Class 0)")
plt.scatter(X[neigh_idx][y[neigh_idx] == 1, 0],
X[neigh_idx][y[neigh_idx] == 1, 1],
s=20, label="Neighbours (Class 1)")
# highlight x0
plt.scatter(x0[0], x0[1],
s=60, color="red", label=r"$x_0$")
plt.xlabel(r"$x_1$")
plt.ylabel(r"$x_2$")
plt.legend()
plt.title(r"Local point $x_0$ and the 50 nearest neighbours")
plt.show()
We can now add synthetic points to further explore the model classifications. We sample our synthetic points from the region around \(x_0\) with the half the variance as the \(J-\) nearest neighbours.
var = np.var(X[neigh_idx])
synthetic_samples = np.random.normal(loc=x0, scale=np.sqrt(var)/2, size=(10_000, len(x0)))
Put our synthetic samples and our real samples into the model to see how they are classified.
real_inputs = X[neigh_idx]
real_preds = model_func(torch.from_numpy(real_inputs).float().to(device))
synthetic_preds = model_func(torch.from_numpy(synthetic_samples).float().to(device))
real_preds.shape
(50,)
plt.scatter(synthetic_samples[synthetic_preds == 0, 0],
synthetic_samples[synthetic_preds == 0, 1],
s=50, facecolor="blue", label="Synthetic samples (class 0)")
plt.scatter(synthetic_samples[synthetic_preds == 1, 0],
synthetic_samples[synthetic_preds == 1, 1],
s=50, facecolor="orange", label="Synthetic samples (class 1)")
plt.scatter(real_inputs[real_preds == 0, 0],
real_inputs[real_preds == 0, 1],
s=50, facecolor="blue", edgecolor="black",
linewidth=1.5, label="Real neighbours (class 0)"
)
# class 1 neighbours
plt.scatter(real_inputs[real_preds == 1, 0],
real_inputs[real_preds == 1, 1],
s=50, facecolor="orange", edgecolor="black",
linewidth=1.5, label="Real neighbours (class 1)"
)
# 4 — highlight x0
plt.scatter(x0[0], x0[1],
s=80, color="red", edgecolor="black", linewidth=1.2,
label=r"$x_0$")
plt.legend()
plt.xlabel(r"$x_1$")
plt.ylabel(r"$x_2$")
plt.title(r"Real and synthetic points in neighbourhood of $x_0$")
Text(0.5, 1.0, 'Real and synthetic points in neighbourhood of $x_0$')
This is how the model classfies the points in this neighbourhood. Now we can approximate the behaviour of this model using a symbolic expression.
from symtorch import SymbolicModel
Detected IPython. Loading juliacall extension. See https://juliapy.github.io/PythonCall.jl/stable/compat/#IPython
Set up the SLIME model.
slime_model = SymbolicModel(model_func, block_name = "func")
slime_params = {"x": x0,
"J_nn": 50,
"num_synthetic": 10_000,
"nn_metric": "euclidean"
}
pysr_params={
'niterations': 1000,
'binary_operators': ['+', '*', '-', '/', 'cond'],
'unary_operators': ['sin', 'cos', 'exp', 'log', 'sqrt'],
'complexity_of_operators': {'sin': 3, 'cos': 3, 'exp': 3, 'log': 3, 'sqrt': 2, 'cond': 2},
'constraints':{'exp': 2, 'log':2, 'sin':2, 'cos':2},
'verbosity': 0
}
slime_model.distill(X, SLIME= True, slime_params= slime_params, sr_params = pysr_params)
🔍 SLIME mode: Using 10050 points (50 real + 10000 synthetic)
Point of interest: [8.23135402 3.41537682]
🛠️ Running SR on output dimension 0 of 0
/Users/liz/PhD/SymTorch_project/symtorch_venv/lib/python3.11/site-packages/pysr/sr.py:2811: UserWarning: Note: it looks like you are running in Jupyter. The progress bar will be turned off.
warnings.warn(
/Users/liz/PhD/SymTorch_project/symtorch_venv/lib/python3.11/site-packages/pysr/sr.py:2265: UserWarning: Note: you are running with more than 10,000 datapoints. You should consider turning on batching (https://ai.damtp.cam.ac.uk/pysr/options/#batching). You should also reconsider if you need that many datapoints. Unless you have a large amount of noise (in which case you should smooth your dataset first), generally < 10,000 datapoints is enough to find a functional form with symbolic regression. More datapoints will lower the search speed.
warnings.warn(
💡Best equation for output 0 found to be cond((((cos(x1) * 0.04937385) + 0.5404889) * x0) + -2.906935, 1.0).
❤️ SR on func complete.
{0: PySRRegressor.equations_ = [
pick score equation \
0 0.000000 0.9994247
1 0.013140 (-0.77621204 / x0) + 1.0939363
2 9.639104 cond(x0 + -5.8403883, 1.0)
3 0.000001 cond(x0 + -5.84091, 0.99998724) - -1.2700881e-5
4 0.063348 cond(((x1 * 0.011053928) + x0) + -5.849988, 1.0)
5 0.061500 cond((((x1 * 0.04135067) + x0) * 1.0558497) + ...
6 0.000024 cond(((x0 * 1.0541468) + (x1 * 0.045976285)) +...
7 >>>> inf cond((((cos(x1) * 0.04937385) + 0.5404889) * x...
loss complexity
0 5.769923e-04 1
1 5.474488e-04 5
2 3.565610e-08 6
3 3.565600e-08 8
4 3.141298e-08 10
5 2.777738e-08 12
6 2.777605e-08 14
7 0.000000e+00 15
]}
This is the equation that the symbolic regression found:
slime_model.show_symbolic_expression(SLIME=True)
➡️ Slime symbolic expressions for output dimension 0:
complexity loss \
0 1 5.769923e-04
1 5 5.474488e-04
2 6 3.565610e-08
3 8 3.565600e-08
4 10 3.141298e-08
5 12 2.777738e-08
6 14 2.777605e-08
7 15 0.000000e+00
equation score \
0 0.9994247 0.000000
1 (-0.77621204 / x0) + 1.0939363 0.013140
2 cond(x0 + -5.8403883, 1.0) 9.639104
3 cond(x0 + -5.84091, 0.99998724) - -1.2700881e-5 0.000001
4 cond(((x1 * 0.011053928) + x0) + -5.849988, 1.0) 0.063348
5 cond((((x1 * 0.04135067) + x0) * 1.0558497) + ... 0.061500
6 cond(((x0 * 1.0541468) + (x1 * 0.045976285)) +... 0.000024
7 cond((((cos(x1) * 0.04937385) + 0.5404889) * x... inf
sympy_format \
0 0.999424700000000
1 1.0939363 - 0.77621204/x0
2 Piecewise((1.0, x0 - 5.8403883 > 0), (0.0, True))
3 Piecewise((0.99998724, x0 - 5.84091 > 0), (0.0...
4 Piecewise((1.0, x0 + x1*0.011053928 - 5.849988...
5 Piecewise((1.0, (x0 + x1*0.04135067)*1.0558497...
6 Piecewise((0.9999528, x0*1.0541468 + x1*0.0459...
7 Piecewise((1.0, x0*(cos(x1)*0.04937385 + 0.540...
lambda_format
0 PySRFunction(X=>0.999424700000000)
1 PySRFunction(X=>1.0939363 - 0.77621204/x0)
2 PySRFunction(X=>Piecewise((1.0, x0 - 5.8403883...
3 PySRFunction(X=>Piecewise((0.99998724, x0 - 5....
4 PySRFunction(X=>Piecewise((1.0, x0 + x1*0.0110...
5 PySRFunction(X=>Piecewise((1.0, (x0 + x1*0.041...
6 PySRFunction(X=>Piecewise((0.9999528, x0*1.054...
7 PySRFunction(X=>Piecewise((1.0, x0*(cos(x1)*0....
🏆 Best: cond((((cos(x1) * 0.04937385) + 0.5404889) * x0) + -2.906935, 1.0) (loss: 0.000000e+00)
which is some sort of conditional decision boundary.
Now let’s see if our SLIME model has approximated the local behaviour well…
slime_predictions=slime_model(synthetic_samples)
slime_labels = (slime_predictions >= 0.5).astype(int)
plt.scatter(
synthetic_samples[slime_labels == 0, 0],
synthetic_samples[slime_labels == 0, 1],
s=50, facecolor="blue", label="Synthetic samples (class 0)"
)
plt.scatter(
synthetic_samples[slime_labels == 1, 0],
synthetic_samples[slime_labels == 1, 1],
s=50, facecolor="orange", label="Synthetic samples (class 1)"
)
plt.legend()
plt.xlabel(r"$x_1$")
plt.ylabel(r"$x_2$")
plt.title(r"SLIME surrogate model predictions of synthetic data in $x_0$ neighbourhood")
Text(0.5, 1.0, 'SLIME surrogate model predictions of synthetic data in $x_0$ neighbourhood')
slime_model.switch_to_symbolic(SLIME=True)
✅ Successfully switched func to SLIME symbolic equations for all 1 dimensions:
Dimension 0: cond((((cos(x1) * 0.04937385) + 0.5404889) * x0) + -2.906935, 1.0)
Variables: ['x0', 'x1']
🎯 All 1 output dimensions now using SLIME symbolic equations.
func = slime_model.get_symbolic_function(SLIME=True)
It has done a pretty good job!
Let’s visualise the SLIME decision boundary on the true dataset.
# make a grid covering data
x0_min, x0_max = X[:,0].min() - 1, X[:,0].max() + 1
x1_min, x1_max = X[:,1].min() - 1, X[:,1].max() + 1
xx, yy = np.meshgrid(
np.linspace(x0_min, x0_max, 600),
np.linspace(x1_min, x1_max, 600)
)
# evaluate model
grid = np.column_stack([xx.ravel(), yy.ravel()]) # shape (600*600, 2)
Z = slime_model(grid)
Z = np.asarray(Z).reshape(-1) # ensure 1D
Z = Z.reshape(xx.shape)
# plot ONLY the boundary line
plt.contour(xx, yy, Z, levels=[0.5], colors='black', linewidths=2)
#plot the original data
plt.scatter(synthetic_samples[synthetic_preds == 0, 0],
synthetic_samples[synthetic_preds == 0, 1],
s=50, facecolor="blue", label="Synthetic samples (class 0)")
plt.scatter(synthetic_samples[synthetic_preds == 1, 0],
synthetic_samples[synthetic_preds == 1, 1],
s=50, facecolor="orange", label="Synthetic samples (class 1)")
plt.scatter(real_inputs[real_preds == 0, 0],
real_inputs[real_preds == 0, 1],
s=50, facecolor="blue", edgecolor="black",
linewidth=1.5, label="Real neighbours (class 0)")
plt.scatter(real_inputs[real_preds == 1, 0],
real_inputs[real_preds == 1, 1],
s=50, facecolor="orange", edgecolor="black",
linewidth=1.5, label="Real neighbours (class 1)")
plt.scatter(x0[0], x0[1],
s=80, color="red", edgecolor="black", linewidth=1.2,
label=r"$x_0$")
plt.xlim(2,15)
plt.ylim(-2,10)
plt.legend()
plt.xlabel(r"$x_1$")
plt.ylabel(r"$x_2$")
plt.title("Decision Boundary of Local Symbolic Surrogate Model")
Text(0.5, 1.0, 'Decision Boundary of Local Symbolic Surrogate Model')
