Additive Gaussian Processes

In plain words

Simple regression models add up one-variable effects. They are easy to read, but they miss how inputs work together. Flexible models let every input affect the output at once. They fit well, but they can be unreliable on combinations of inputs they have not seen.

This paper's model sits between the two. Its prior over functions is a sum of pieces. Each piece depends on a small group of inputs: single inputs, pairs, triples, and so on, up to all inputs together. The model learns how much weight each level of interaction gets, so the fit shows which kind of structure the data has. Summing every group naively would be very expensive, but a known algebra identity reduces the work to a cost that grows with the square of the number of inputs. The authors report results that are the best or statistically tied on most of the datasets they test.

TL;DR

• The kernel is a sum over interaction orders n = 1..D of an order-n term, each scaled by its own variance σ_n² §3.1 eq. 1–eq. 4.
• Base kernels are one-dimensional, one per input; with squared-exponential bases the top order is the usual SE kernel eq. 5.
• All orders are summed in O(D²) with Newton–Girard identities §3.3 eq. 6.
• Order variances show which interaction orders explain the target §3.2 tab. 1.
• Results are author-reported on 5 regression and 6 classification datasets tab. 2–tab. 5.
• The official code is GPL-2.0: cite it by path and line only.

Problem

• Generalized additive models (GAMs) sum one-dimensional functions. They are easy to fit and interpret, but they cannot express interactions between inputs §1.
• Squared-exponential (SE) GPs let all inputs act at once. They are flexible but local, so they generalize poorly to new combinations of inputs, and they need data at all relevant input combinations §1 §4.3.
• Earlier methods that allow interactions either select a hull of interaction terms by cross-validation (hierarchical kernel learning, HKL) §4.1, restrict SS-ANOVA to first and second order in practice §4.2, or use one interaction order only (support-vector ANOVA) §4.2.
• Kernel choice has a large effect on modeling quality for small to medium datasets §2.

Method

The model is a Gaussian process whose covariance is an additive kernel built from one-dimensional base kernels. The only design choice is one base kernel per input dimension. Their hyperparameters, such as length-scales, are learned by maximizing the marginal likelihood. The model adds one variance per interaction order. Everything else in GP inference is unchanged §3.1 §3.4.

1. Choose a one-dimensional base kernel k_d(x_d, x'_d) for each input dimension d = 1..D. This is the only design choice §3.1.
2. Define the order-n additive kernel as the variance σ_n² times the sum over all n-subsets of dimensions of the product of their base kernels. Each order has its own σ_n² §3.1 eq. 1–eq. 3.
3. The order-D term is one product over all dimensions eq. 4.
4. With SE base kernels, the order-D term is the multivariate SE (Gaussian) kernel eq. 5.
5. The full kernel is the sum of orders 1 to D §3.1. The number of terms in order n is C(D, n) §3.1.
6. Efficient evaluation: set z_d = k_d(x_d, x'_d). Order n is the n-th elementary symmetric polynomial e_n(z_1..z_D). The Newton–Girard recursion with e_0 = 1 computes all orders in O(D²) §3.3 eq. 6.
7. Derivatives with respect to a base kernel: ∂k_add_n / ∂z_j = e_{n-1} over all dimensions except j eq. 7 §3.3. This makes the base-kernel hyperparameters differentiable for the marginal likelihood §3.3.
8. Gram-matrix cost is O(N²DR) with R the maximum interaction order. The authors note that the first few orders usually matter most, so R can be lowered with little loss §3.4.
9. Fit: maximize the training marginal likelihood over the base-kernel hyperparameters and the order variances σ_1..σ_D, using a constant mean §5.2.

k_{add_n}(\mathbf{x}, \mathbf{x}') = \sigma_n^2 \sum_{1 \le i_1 < \dots < i_n \le D} \prod_{d=1}^{n} k_{i_d}(x_{i_d}, x'_{i_d}) eq. 1–eq. 3

e_n(z_1..z_D) = \frac{1}{n} \sum_{k=1}^{n} (-1)^{k-1} \, e_{n-k}(z_1..z_D) \, s_k(z_1..z_D), \qquad s_k = \sum_{d=1}^{D} z_d^{\,k}, \quad e_0 = 1 eq. 6

Kernel evaluation over the orders, in our own notation:

# Additive kernel over D inputs, orders 1..R  [§3.1] [§3.3] [eq. 1]-[eq. 3]
# Input:  X1 (N1, D), X2 (N2, D); lengthscales (D,); order_var (R,)
# Output: K (N1, N2), the Gram matrix

def additive_kernel(X1, X2, lengthscales, order_var):
    # One-dimensional SE base kernel per input dimension  [§3.1]
    z = [
        se_kernel(X1[:, d], X2[:, d], lengthscales[d])
        for d in range(D)
    ]
    # All elementary symmetric polynomials e_0..e_R of z  [§3.3] [eq. 6]
    e = elementary_symmetric(z, R)
    # Order n is weighted by its own variance sigma_n^2  [eq. 1]-[eq. 3]
    K = zeros(N1, N2)
    for n in range(1, R + 1):
        K = K + order_var[n - 1] * e[n]
    return K

Newton–Girard recursion, in our own notation:

# Elementary symmetric polynomials by Newton-Girard  [§3.3] [eq. 6]
# e_0 = 1; for n = 1..R:
#   e_n = (1/n) * sum over k = 1..n of (-1)^(k-1) * e_(n-k) * s_k
# s_k = sum over dimensions d of z_d^k  (power sums)

def elementary_symmetric(z, R):
    # Power sums s_1..s_R over the D base kernels  [eq. 6]
    power_sums = [sum(zd ** k for zd in z) for k in range(1, R + 1)]
    e = [ones_like(z[0])]  # e_0 = 1  [eq. 6]
    for n in range(1, R + 1):
        total = zeros_like(z[0])
        for k in range(1, n + 1):
            total = total + (-1) ** (k - 1) * e[n - k] * power_sums[k - 1]
        e.append(total / n)  # e_n  [eq. 6]
    return e

Training procedure, in our own notation (initialisation details are in the hyperparameters table and are not stated in the paper):

# Fit by maximizing the log marginal likelihood  [§5.2]
# Input:  X (N, D), y (N,); R = min(D, 10)  [§5.2]
# Output: fitted lengthscales, order variances, noise, mean

def fit(X, y):
    # Input standardisation used in the code, not stated in the paper
    X = (X - X.mean(axis=0)) / X.std(axis=0)  # [code: run_one_fold.m#L40-L41]
    # Initial values: see Hyperparameters table  [code: gp_add_general_lo.m#L15-L43]
    theta = initial_hyperparameters(D, R)
    # L-BFGS for 500 iterations on the negative log marginal likelihood  [§5.2]
    theta = lbfgs(neg_log_marginal_likelihood, theta, iterations=500)
    return theta

Hyperparameters

Results

Author-reported. Setup: mean across 10 train-test splits §5.3 tab. 2–tab. 5. Regression: five datasets (bach, concrete, pumadyn-8nh, servo, housing), with linear regression, GP-GAM (first order only), HKL, SE-GP with ARD, and the additive GP §5.2. Classification: six datasets (breast, pima, sonar, ionosphere, liver, heart), with logistic regression, GP-GAM, HKL, SE-GP and the additive GP §5.2. HKL has no noise model, so it is left out of the negative-log-likelihood tables §5.3. In the paper, bold means best or not significantly different from the best (paired t-test) §5.3.

Order-variance shares (percent of variance per order, max order 10) tab. 1

Regression, mean squared error (lower is better) tab. 2

Regression, negative log likelihood (lower is better; HKL not reported) tab. 3

Classification, percent error (lower is better) tab. 4

Classification, negative log likelihood (lower is better) tab. 5

Per §5.3, no competitor significantly beats the additive model, and on some datasets it beats all of them significantly §5.3. The tables agree on regression; on heart the additive GP is not in the bold set for error tab. 4. The authors attribute slightly worse cases to overfitting or to a local maximum in hyperparameter search §5.3.

Figures, described in words: fig. 1 draws from a first-order additive kernel versus a product kernel. fig. 2 shows low-order functions recovered on the concrete dataset. fig. 3 compares HKL hulls, GP-GAM, SE-GP and the additive kernel. fig. 4 shows kernel isocontours in three dimensions. fig. 5 is a synthetic sum of two axis-aligned sine functions, trained on an L-shaped region and tested on a peak away from it; the additive GP recovers both sine functions and most of the variance is first order §5.1.

Paper vs code

Official repository: duvenaud/additive-gps (owner is David Duvenaud, the first author; the README cites this arXiv id and links his homepage), commit 534edd1b04ad670dd713244c1796cf80f093b091, license GPL-2.0. Per the rules for copyleft code, no code is quoted. Everything below is cited by path and line.

• Optimizer budget. Paper: L-BFGS for 500 iterations §5.2. Code: minimize is called with −500, and the header says a negative count means function evaluations gp_add_general_lo.m gp_add_general_lo.m minimize.m. So the code budgets 500 function evaluations, not 500 iterations.
• Restarts. Paper: five random restarts §5.2. Code: the compiler's default seeds are 0–4 and the best run by training marginal likelihood is kept compile_all_results.m compile_results.m. But the driver calls compile_all_results(K, seed, outdir) with a scalar seed of 0, so the default is not used there run_all_experiments.m run_all_experiments.m. The card cannot tell which call produced tables 2–5.
• Folds. Paper: mean across 10 train-test splits §5.3. Code: K = 10, but the permutation is the identity, so folds are contiguous blocks in file row order, not random splits run_all_experiments.m run_one_fold.m gen_kfolds.m. The seed controls only the random initialisation.
• Dataset files. Table 2 names "pumadyn-8nh" and "concrete". The code's regression list names r_pumadyn512.mat, r_concrete_500, bach_synth_r_200.mat, r_servo.mat and r_housing.mat define_datasets_and_methods.m. The card cannot confirm these are the paper's exact files; the source does not say how the r_ files were made.
• Standardisation. Not in the paper. Code standardises every input column; regression targets are also standardised, classification labels are not run_one_fold.m.
• Base kernel. Paper: base-kernel hyperparameters are learned, the parameterisation is not specified §3.1. Code: one length-scale per dimension and no per-dimension variance covSEiso_length.m.
• Order-variance parameterisation. Code stores log σ_n and covADD computes the order variance as exp(2·h) covADD.m. The initialisation writes log(v) into h gp_add_general_lo.m, so the effective starting variance appears to be v², not v. This was read from code, not run.
• Hyper-prior penalty. Code adds Σ(h_i/10)^10 to the negative log likelihood penalized_gp.m. Not in the paper.
• Consistent. Max order 10 §5.2 common_gp_parameters.m; constant mean §5.2 gp_add_general_lo.m; EP for classification §5.2 gp_add_class_lo.m; GP-GAM as the max-order-1 variant gp_add_lo_1.m.

Pitfalls

• Initialisation is not in the paper. Length-scales, order variances and noise are random or fixed in code (see Hyperparameters). Results depend on seed; the code keeps the run with the highest training marginal likelihood. The paper does not say how restarts are chosen.
• Stored log-variance is squared on use. Follow the code only if you reproduce its scale: a faithful implementation of the paper's σ_n² would not square the initial value.
• Order-variance scale is combinatorial. The order-n term is a sum of C(D, n) products, so raw magnitudes differ by orders of magnitude. The code divides the initial order variances by the order-n kernel value at zero distance gp_add_general_lo.m. The paper does not discuss this.
• Numerical stability is not discussed. The Newton–Girard recursion alternates signs eq. 6 and can lose precision at large D or large R. This card does not test it.
• Gradient cost. For each base-kernel hyperparameter, the code recomputes a leave-one-dimension-out polynomial, so the gradient cost grows with D on top of the kernel cost covADD.m. The paper states only the forward cost §3.4.
• Interpretation depends on the fit. Order variances describe the fitted model, not the data alone. The authors present them as a way to find which orders matter §1 §3.2 tab. 1.
• Axis-aligned by construction. Interactions are between input coordinates, so rotated structure is not captured without rotating the inputs. The compiled paper does not discuss this limit; the card infers it from the definition §3.1.
• Hyperparameters are point estimates. The authors expect integrating them out to help §5.3.
• Bach dataset. The compiled paper does not describe the "bach" dataset, so this card does not describe it either. Some text in the source comes after \end{document} and is not part of the paper.
• Author-reported only. Tables are 10-fold means with paired t-tests as stated in the paper. No replication was run.

Lineage

• Builds on: GPs and squared-exponential kernels §2. Hierarchical kernel learning (HKL, Bach), which this paper compares against in its experiments §4.1. Plate's additive GP with first and D-th order terms only §4. Support-vector ANOVA kernels §4.2 and smoothing-splines ANOVA §4.2.
• Compared against: HKL, GP-GAM, SE-GP with ARD, linear regression, logistic regression §5.2 tab. 2–tab. 5.
• Superseded by: not known from content already in this wiki.