Lecture 14, Notebook 07: Active subspaces — a nonlinear 10D function - Deep Learning for Solving and Estimating Dynamic Models in Economics and Finance

Course: Deep Learning for Solving and Estimating Dynamic Models in Economics and Finance
Script reference: §9.5 (Active subspaces — nonlinear 10D extension where linear AS needs $d \\geq 2$ )
Notebook role: extension
Author: Simon Scheidegger

RUN_MODE = "smoke"  # one of: "smoke", "teaching", "production"
SEED = 0

The previous notebooks showed functions where a single active dimension ( $d=1$ ) suffices. Here we introduce a nonlinear 10D test function that involves a product term $x_2 \cdot x_3$ , making the gradient depend on the input location. As a result, a 1D active subspace is no longer sufficient -- we need $d = 2$ or $d = 3$ active dimensions to achieve good surrogate accuracy. This notebook compares ASGP surrogates with $d \in \{1, 2, 3\}$ active dimensions against a full 10D GP.

Reference: Scheidegger & Bilionis (2019), Machine Learning for High-Dimensional Dynamic Stochastic Economies, Journal of Computational Science 33, 68--82.

Extends: The examples in Section 3 of the paper with a nonlinear variant that requires multiple active dimensions.

import numpy as np
import matplotlib.pyplot as plt
from scipy import linalg
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C

plt.rcParams['font.size'] = 13

# Hyperparameter budget, dispatched on RUN_MODE (see the run-mode cell above).
if RUN_MODE == "smoke":
    N_GRAD, N_RESTARTS = 200, 0
elif RUN_MODE == "teaching":
    N_GRAD, N_RESTARTS = 500, 3
elif RUN_MODE == "production":
    N_GRAD, N_RESTARTS = 2_000, 10
else:
    raise ValueError(f"Unknown RUN_MODE={RUN_MODE!r}")

1. Nonlinear Test Function¶

We define a nonlinear variant of the 10D exponential function by multiplying with $x_2 \cdot x_3$ :

f(x) = x_2 \cdot x_3 \cdot \exp\!\bigl(0.01\,x_1 + 0.7\,x_2 + 0.02\,x_3 + \ldots + 0.1\,x_{10}\bigr)

(1)

The product $x_2\,x_3$ introduces cross-terms in the gradient: $\partial f / \partial x_2$ depends on $x_3$ and vice versa. This means the gradient direction is no longer constant across input space, so a single linear projection cannot capture all the variation. We also define the baseline exponential test_old (used internally for gradient computation).

def test_old(x):
    return np.exp(0.01*x[0] + 0.7*x[1] + 0.02*x[2] + 0.03*x[3] + 0.04*x[4] + 
                  0.05*x[5] + 0.06*x[6] + 0.08*x[7] + 0.09*x[8] + 0.1*x[9])


def test_function(x):
    return x[1]*x[2] * np.exp(0.01*x[0] + 0.7*x[1] + 0.02*x[2] + 0.03*x[3] + 0.04*x[4] + 
                  0.05*x[5] + 0.06*x[6] + 0.08*x[7] + 0.09*x[8] + 0.1*x[9])

def dtest_function(x):
    _test_old = test_old(x)
    val = np.atleast_1d(x[1] * x[2] * _test_old)
    coefs = np.array([0.01, 0.7, 0.02, 0.03, 0.04, 0.05, 0.06, 0.08, 0.09, 0.1])
    out = val[:, None] * coefs[None, :]
    out[:, 1] += x[2] * _test_old
    out[:, 2] += x[1] * _test_old
    return out

# Random points on Omega = [-1,1]^D
def randOmega(N, D):
    "random points on \Omega = [-1,1]^D"
    return 2 * (np.random.rand(N, D) - 0.5)

2. Active Subspace Construction¶

For each training-set size $N \in \{10, 30, 100, 250, 500, 1000\}$ , we:

Sample $N$ points in $[-1,1]^{10}$ and evaluate $f$ and $\nabla f$ .
Build $C_N = \frac{1}{N} G^\top G$ and eigen-decompose it.
For each active subspace dimension $d \in \{1, 2, 3\}$ , project onto the top- $d$ eigenvectors $W_d \in \mathbb{R}^{10 \times d}$ and fit an ASGP on the $d$ -dimensional projected inputs.
Also fit a standard GP on the full 10D inputs for comparison.
Report the maximum absolute error on 1,000 held-out test points.

Note on GP hyperparameter optimization. Below we set n_restarts_optimizer=0 on the ASGP fits to keep the runtime of this sweep modest in the classroom. For production use, set n_restarts_optimizer >= 3 to avoid local minima in length-scale optimization; we use 0 here only for runtime.

def test_case():
    np.random.seed(SEED)
    Nvals = np.array([10, 30, 100, 250, 500, 1000])
    num_Ns = len(Nvals)
    max_errs = np.inf * np.ones((num_Ns, 3))
    max_errs_gp = np.inf * np.ones(num_Ns)
    
    # construct test points
    X_test = randOmega(1000, 10)
    f_test = test_function(X_test.T)

    for i, N in enumerate(Nvals):
        # training points
        X = randOmega(int(N), 10)
        V = test_function(X.T)
        G = dtest_function(X.T)
    
        CN = (G.T @ G) / N
        vals, vecs = linalg.eigh(CN)
        for d in range(3, 0, -1):

            # find active subspace of dimension d
            W = vecs[:, -d:]
            Y = X @ W

            # fit GP on active subspace
            gp_as = GaussianProcessRegressor(RBF(), n_restarts_optimizer=N_RESTARTS)
            gp_as.fit(Y.reshape(N, d), V)
            m_tilde = gp_as.predict((X_test @ W).reshape(1000, d))
            max_errs[i, d-1] = np.max(np.abs(f_test - m_tilde))
        
        gp = GaussianProcessRegressor(RBF())
        gp.fit(X, V)
        m_tilde = gp.predict(X_test)
        max_errs_gp[i] = np.max(np.abs(f_test - m_tilde))
    
    return vals, Nvals, max_errs, max_errs_gp

vals, Nvals, max_errs, max_errs_gp = test_case()

3. Visualization¶

Left panel -- Eigenvalue spectrum: Unlike the purely linear-exponent case in Notebook 06, the eigenvalue spectrum here shows a more gradual decay. The top eigenvalue still dominates, but the second and third eigenvalues are non-negligible, reflecting the additional directions introduced by the $x_2 \cdot x_3$ product term.

Right panel -- Max absolute error vs. training size $N$ :

1D AS (single active dimension): Performs reasonably but plateaus because it cannot capture the interaction between $x_2$ and $x_3$ .
2D AS: A substantial improvement -- the second active dimension captures the cross-term structure.
3D AS: Further improvement, achieving accuracy close to or better than the full 10D GP while operating in a much lower-dimensional space.
Full GP (10D): Suffers from the curse of dimensionality at small $N$ but improves with more data.

The key lesson: when $d=1$ is not enough, the eigenvalue spectrum tells you how many active dimensions to use.

fig, ax = plt.subplots(1, 2, figsize=(9, 5))
x = np.arange(1, 11)
ax[0].semilogy(x, vals, ".", ms=15)
ax[0].set_xlabel("Eigenvalues")
ax[0].set_xticks(x)
ax[0].set_ylabel(r"$\lambda$")

ax[1].semilogy(Nvals, max_errs, ".-", ms=15)
ax[1].semilogy(Nvals, max_errs_gp, ".-", ms=15)
ax[1].legend([str(i) + "d AS" for i in range(1, 4)] + ["Full GP"])
ax[1].set_xlabel("# of points")

fig.tight_layout()
plt.show()

4. Take-Aways¶

The cross-term $x_2 x_3$ mixes two input directions so that the gradient direction varies across input space. The gradient outer product $C_N$ therefore picks up more than one nonzero eigendirection.
Adding a second active direction ( $d = 2$ ) yields a substantial drop in surrogate error; a third direction buys a bit more. The elbow of the eigenvalue spectrum aligns with the elbow of the error curve: the spectrum is a predictive diagnostic for how many active dimensions are worth keeping.
The full 10D GP eventually catches up as $N$ grows, but for small- $N$ (the regime that matters for expensive simulators) the ASGP with the right $d$ is markedly more data-efficient.

When even $d = 2$ is too many, the active manifold may be curved rather than linear: two linear features combined through a nonlinear aggregator. Notebook 09_Deep_Active_Subspace_Ridge shows a constructed target where linear AS requires $d = 2$ but a deep encoder collapses the same problem to $d = 1$ .

Takeaway¶

When the gradient direction varies across the input space (the $x_2 x_3$ product term), linear AS needs $d = 2$ , not $d = 1$ . The deep-AS recipe of NB 09–10 generalises further to curved manifolds.

# Smoke-mode validation: the active-subspace sweep produced finite eigenvalues and finite errors.
if RUN_MODE == 'smoke':
    assert np.isfinite(vals).all(), 'Active-subspace eigenvalues are not all finite.'
    assert np.isfinite(max_errs).all() and np.isfinite(max_errs_gp).all(), 'AS / GP errors are not all finite.'