LIME from Scratch in C++

July 2026 · 18 min read

This post implements the core of LIME for numeric tabular rows. It samples perturbations around a row, queries the black-box predictor, weights each perturbation by distance, then solves weighted linear regression. The resulting coefficients explain the prediction locally.

Weighted Regression Helper

using Vector = vector<double>;
using Matrix = vector<Vector>;
using Predictor = function<double(const Vector&)>;

Vector solve_linear(Matrix A, Vector b) {
    int n = (int)b.size();
    for (int col = 0; col < n; ++col) {
        int piv = col;
        for (int r = col + 1; r < n; ++r)
            if (fabs(A[r][col]) > fabs(A[piv][col])) piv = r;
        swap(A[col], A[piv]); swap(b[col], b[piv]);
        double div = A[col][col];
        for (int c = col; c < n; ++c) A[col][c] /= div;
        b[col] /= div;
        for (int r = 0; r < n; ++r) if (r != col) {
            double m = A[r][col];
            for (int c = col; c < n; ++c) A[r][c] -= m * A[col][c];
            b[r] -= m * b[col];
        }
    }
    return b;
}

LIME Core

Vector lime_explain(const Predictor& pred,
                    const Vector& x0,
                    int samples,
                    double noiseStd,
                    double kernelWidth,
                    mt19937& rng) {
    int d = (int)x0.size();
    normal_distribution<double> noise(0.0, noiseStd);

    Matrix xtx(d + 1, Vector(d + 1, 0.0));
    Vector xty(d + 1, 0.0);

    for (int s = 0; s < samples; ++s) {
        Vector z = x0;
        double dist2 = 0.0;
        for (int j = 0; j < d; ++j) {
            double eps = noise(rng);
            z[j] += eps;
            dist2 += eps * eps;
        }
        double w = exp(-dist2 / (kernelWidth * kernelWidth));
        double y = pred(z);

        Vector row(d + 1, 1.0);              // intercept + features
        for (int j = 0; j < d; ++j) row[j + 1] = z[j] - x0[j];

        for (int a = 0; a < d + 1; ++a) {
            xty[a] += w * row[a] * y;
            for (int b = 0; b < d + 1; ++b)
                xtx[a][b] += w * row[a] * row[b];
        }
    }

    for (int j = 0; j < d + 1; ++j) xtx[j][j] += 1e-6; // ridge stabilizer
    Vector beta = solve_linear(xtx, xty);
    return Vector(beta.begin() + 1, beta.end());        // feature coefficients
}