Kernel SHAP from Scratch in C++

July 2026 · 20 min read

Kernel SHAP estimates Shapley values by fitting a weighted linear model on coalition masks. The black-box model is queried on rows where some features come from the row being explained and absent features come from a background row.

Coalition Helpers

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

double choose_weight(int m, int s) {
    if (s == 0 || s == m) return 1000.0;   // anchor empty and full coalitions
    double comb = 1.0;
    for (int i = 1; i <= s; ++i) comb = comb * (m - s + i) / i;
    return (m - 1.0) / (comb * s * (m - s));
}

Vector masked_row(const Vector& x0, const Vector& bg, const vector<int>& mask) {
    Vector z = x0;
    for (size_t j = 0; j < z.size(); ++j)
        if (!mask[j]) z[j] = bg[j];
    return z;
}

Kernel SHAP Core

The same solve_linear helper from the LIME post solves the weighted normal equations. The design matrix is now binary masks rather than numeric perturbations.

Vector kernel_shap(const Predictor& pred,
                   const Vector& x0,
                   const Matrix& background) {
    int m = (int)x0.size();
    int cols = m + 1; // intercept + one coefficient per feature
    Matrix xtx(cols, Vector(cols, 0.0));
    Vector xty(cols, 0.0);

    for (int maskBits = 0; maskBits < (1 << m); ++maskBits) {
        vector<int> mask(m, 0);
        int s = 0;
        for (int j = 0; j < m; ++j) {
            mask[j] = (maskBits >> j) & 1;
            s += mask[j];
        }
        double w = choose_weight(m, s);

        for (const auto& bg : background) {
            Vector z = masked_row(x0, bg, mask);
            double y = pred(z);
            Vector row(cols, 1.0);
            for (int j = 0; j < m; ++j) row[j + 1] = mask[j];

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

    for (int j = 0; j < cols; ++j) xtx[j][j] += 1e-8;
    Vector beta = solve_linear(xtx, xty);
    return Vector(beta.begin() + 1, beta.end());
}