Cost-Aware Active Learning with Directional Derivatives

This chapter describes JetGP’s active-learning framework for adaptively choosing where to evaluate an expensive function and along which directions to evaluate its derivatives. The framework’s distinguishing feature is a single unified, cost-aware acquisition rule that compares function and directional-derivative candidates against each other in normalised information-per-cost units, rather than treating direction selection as a separate post-processing step after a function point has been chosen.

Notation follows the rest of the theory manual: \(f:\mathcal{X}\to\mathbb{R}\) is the expensive scalar response over a bounded input domain \(\mathcal{X}\subset\mathbb{R}^d\), \(p(\mathbf{x})\) is the input density associated with the task (uniform if unspecified), and the GP model follows the directional-derivative formulation summarised in Directional Derivative-Enhanced Gaussian Processes.

—

Acquisition Quantities

After the current GP has been refit, two families of candidate observations are scored: a single new function evaluation somewhere in \(\mathcal{X}\), and one or more new directional derivatives at the existing training points.

Function candidate

The next candidate function location is obtained by maximising the PDF-weighted log-variance acquisition over the bounded domain:

\[\mathbf{x}_{\mathrm{new}} \;=\; \arg\max_{\mathbf{x}\in\mathcal{X}}\bigl[\; \log \sigma_f^2(\mathbf{x}) \;+\; \log p(\mathbf{x}) \;\bigr],\]

which is equivalent to maximising \(\sigma_f^2(\mathbf{x})\,p(\mathbf{x})\) wherever the density is positive. The log form is used for numerical stability when the posterior variance or density weight is small. When \(p\) is uniform this reduces to plain maximum posterior variance.

The optimisation is performed by Latin-hypercube multistart followed by local bounded L-BFGS-B from each start. A large finite penalty is returned for infeasible or degenerate evaluations, so the finite-difference gradient needed by L-BFGS-B remains well-defined.

Directional-derivative candidates

At each existing training point \(\mathbf{x}_i\), the posterior gradient covariance under the current GP is

\[\mathbf{C}_d(\mathbf{x}_i) \;=\; \mathrm{Cov}\bigl[\nabla f(\mathbf{x}_i)\bigr] \in \mathbb{R}^{d\times d}.\]

For any unit direction \(\mathbf{v}\), the posterior variance of the directional derivative is the Rayleigh quotient

\[\sigma^2_{\partial_{\mathbf{v}} f}(\mathbf{x}_i) \;=\; \mathbf{v}^{\top}\mathbf{C}_d(\mathbf{x}_i)\mathbf{v}.\]

Maximising this over \(\|\mathbf{v}\|=1\) yields the dominant eigenvector of \(\mathbf{C}_d(\mathbf{x}_i)\), with the maximum value being the largest eigenvalue (Rayleigh-Ritz). Subsequent directions, taken over the orthogonal complement of previously chosen directions at the same anchor, are the remaining eigenvectors in descending eigenvalue order.

Each candidate direction is normalised against a prior reference scale,

\[\rho_j \;=\; \frac{\lambda_j\!\bigl(\mathbf{C}_d(\mathbf{x}_i)\bigr)} {\lambda_{\max}\!\bigl(\mathbf{K}_{\nabla}^{\mathrm{prior}}\bigr)},\]

where \(\mathbf{K}_{\nabla}^{\mathrm{prior}} = \mathbf{K}_{\nabla\nabla}(\mathbf{x},\mathbf{x})\) is the gradient covariance implied by the kernel before conditioning on data. The ratio \(\rho_j\in[0,1]\) measures the residual posterior uncertainty in direction \(\mathbf{v}_j\) as a fraction of the largest derivative uncertainty the model can express.

Implementation note (Lanczos top-\(k\))

Full eigendecomposition of \(\mathbf{C}_d(\mathbf{x}_i)\) is avoided in the implementation. The \(k\) leading eigenpairs are obtained by Lanczos iteration with full re-orthogonalisation, where \(k = \min(\texttt{max\_directions},\, d) - q_i\) and \(q_i\) is the number of directional observations already taken at \(\mathbf{x}_i\). This is significantly cheaper for high-dimensional inputs and naturally respects the per-anchor cap max_directions.

—

Unified Cost-Aware Score

The cost-aware framework departs from the classical two-stage rule (pick \(\mathbf{x}\) by MSE, then pick directions at \(\mathbf{x}\)) by ranking function and derivative candidates against each other under a single information-per-cost score.

Let \(c_f\) and \(c_d\) denote the costs of one function evaluation and one directional-derivative observation respectively (both in arbitrary matching units; see the discussion below). The score for the function candidate at \(\mathbf{x}_{\mathrm{new}}\) is

\[s_f \;=\; \frac{\sigma_f^2(\mathbf{x}_{\mathrm{new}})} {\lambda_{\max}\!\bigl(\mathbf{K}_f^{\mathrm{prior}}\bigr)} \,\cdot\, p(\mathbf{x}_{\mathrm{new}}) \,\big/\, c_f \;=\; \rho_f\, p(\mathbf{x}_{\mathrm{new}}) \,/\, c_f,\]

and the score for a directional-derivative candidate at anchor \(\mathbf{x}_i\) with direction \(\mathbf{v}_j\) is

\[s_{d,ij} \;=\; \rho_j(\mathbf{x}_i)\, p(\mathbf{x}_i) \,/\, c_d.\]

Both scores are dimensionless (information ratio per unit cost), so they are commensurable across modalities. The next observation is then chosen by joint argmax,

\[\mathrm{candidate}^{*} \;=\; \arg\max\bigl\{\,s_f,\; \{s_{d,ij}\}_{i,j}\,\bigr\},\]

subject to the affordability constraint \(c(\mathrm{candidate}) \le B - C_{\mathrm{spent}}\) where \(B\) is the total cost budget and \(C_{\mathrm{spent}}\) is the cumulative cost spent so far.

A relative tolerance \(\rho_{\mathrm{tol}}\) discards candidates whose underlying \(\rho\) (before density and cost scaling) falls below the gate; setting \(\rho_{\mathrm{tol}}=0\) keeps every feasible candidate.

Interpretation

• When \(c_d \ll c_f\) (e.g. automatic differentiation or OTI), the derivative scores are upweighted by \(c_f/c_d\), and the policy naturally prefers directional derivatives at existing anchors over fresh function evaluations.

• When \(c_d \gg c_f\) (e.g. central finite differences, which cost \(2c_f\) per directional derivative), function evaluations dominate.

• The density factor \(p(\mathbf{x})\) is what makes the framework test-distribution-aware: candidates at low-density inputs are deprioritised even if their posterior uncertainty is high. Without this factor the policy chases variance at the corners of \(\mathcal{X}\), which may have negligible weight under the task distribution.

—

Hyperparameter Optimisation

The kernel hyperparameters \(\boldsymbol{\psi}\) are obtained by maximising the marginal log-likelihood corresponding to the current mixed observation set:

\[\boldsymbol{\psi}^* \;=\; \arg\max_{\boldsymbol{\psi}}\; \log p(\mathbf{y}^{DD}\mid \boldsymbol{\psi}),\]

where \(\mathbf{y}^{DD}\) is the augmented vector of function values and directional derivatives observed so far, and the likelihood is constructed from the directional-derivative covariance described in Directional Derivative-Enhanced Gaussian Processes. JetGP’s hybrid global+local optimiser is used: a JADE adaptive differential-evolution global search followed by periodic L-BFGS-B refinement of the best population member. Each refit after the first is warm-started from the previously optimised hyperparameters, which is appropriate in sequential design where each model update only modestly perturbs the observation set.

—

Algorithm Summary

The adaptive procedure has two stages: an initial design-and-fit, then a single budget-constrained sequential loop with a unified candidate set. Unlike classical two-stage schemes that pick a function point first and then choose directions at that point, the cost-aware implementation rebuilds one mixed candidate set every iteration and picks the single best candidate across modalities and anchors.

Inputs:  initial DOE size n0; cost budget B; iteration cap n_iter;
         per-modality costs c_f, c_d, (c_d2 if 2nd-order is enabled);
         max_directions cap m; tolerance rho_tol; input density p(x);
         acquire_second_order flag.
Output:  a directional-derivative-enhanced GP surrogate.

1. Initialise.
     a. Sample (or accept) an initial DOE X_0 of n0 points.
     b. Evaluate f(X_0).
     c. Fit the GP by maximising the marginal log-likelihood.

2. Sequential design loop. Set cumulative_cost = 0. While
   cumulative_cost < B and step < n_iter:

     a. Build the unified candidate set:
          i.  One function candidate at x_new = argmax of the wMPV
              log-objective over the input bounds.
          ii. For each existing anchor x_i:
              - extract the top-k eigenpairs of C_d(x_i) via Lanczos,
                where k = min(m, d) - q_i and q_i is the number of
                directional observations already taken at x_i,
              - keep eigenvector v_j iff rho_j > rho_tol,
              - emit a 1st-order derivative candidate (kind=d, order=1,
                cost=c_d) for each surviving (x_i, v_j),
              - if acquire_second_order is enabled, also emit a 2nd-order
                derivative candidate (kind=d, order=2, cost=c_d2) for
                each surviving (x_i, v_j) whose 2nd-order rho clears the
                tolerance.

     b. Score every candidate as s = rho * p(x_eval) / c, where x_eval
        is x_new for function candidates and x_i for derivative
        candidates.

     c. Drop candidates whose cost would exceed the remaining budget
        (B - cumulative_cost).

     d. If no candidate survives, exit the loop (no admissible
        candidate above rho_tol).

     e. Pick the candidate with the largest score; evaluate it; update
        the training data; refit the GP (warm-started from the
        previous hyperparameters); update cumulative_cost.

The loop terminates on the earliest of three conditions: (i) the iteration cap n_iter is reached, (ii) the cost budget is exhausted, or (iii) no affordable candidate clears the tolerance. The max_directions cap m is enforced at Lanczos extraction time — only m - q_i additional directions are even considered at an anchor that already has q_i directional observations — which prevents pathological stacking of derivative observations at a single point. When a function candidate is chosen mid-loop, the new training point immediately becomes a valid anchor for derivative candidates on subsequent iterations.

—

Relationship to Classical Active Learning Criteria

The cost-aware framework is a strict generalisation of the classical maximum-posterior-variance (MPV) and PDF-weighted MPV criteria:

• If \(c_d=\infty\) (derivatives forbidden) and \(p\) is uniform, the function-only score reduces to plain MPV.

• If \(c_d=\infty\) and \(p\not\equiv\mathrm{const.}\), the function-only score reduces to weighted MPV (the pointwise greedy approximation of weighted IMSE).

• If \(c_f=\infty\) (function observations forbidden), the policy becomes a pure-directional-derivative active learner, picking eigenvectors of \(\mathbf{C}_d\) at existing anchors.

• For finite \(c_f, c_d\), the policy adaptively blends both modalities in proportion to their normalised information-per-cost ratio.

—

Limitations and Future Extensions

The acquisition score is local and greedy. Each \(\rho\) measures posterior uncertainty reduction at the anchor point itself, not at the distribution of test points where prediction accuracy will ultimately be evaluated. This is mathematically equivalent to a one-step look-ahead in normalised posterior variance per cost, but only an approximation of the true objective (integrated test-distribution mean-squared error).

The gap between the two is small for well-conditioned input distributions with a clear density gradient, but it can widen for mixed input distributions with flat-density marginals, where the cost-aware framework can over-commit to a small number of high-variance anchor points. A rigorous fix would replace the local \(\rho\,p(\mathbf{x})\) score by the expected reduction in integrated posterior variance over a Monte-Carlo representation of the test distribution (an IMSE / predictive-variance- reduction criterion). This is a planned future extension.

—

References

[1]

A Comparison of Numerical Optimizers in Developing High Dimensional Surrogate Models, volume Volume 2B: 45th Design Automation Conference of International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, 08 2019. URL: https://doi.org/10.1115/DETC2019-97499, arXiv:https://asmedigitalcollection.asme.org/IDETC-CIE/proceedings-pdf/IDETC-CIE2019/59193/V02BT03A037/6452976/v02bt03a037-detc2019-97499.pdf, doi:10.1115/DETC2019-97499.

[2]

Matheron Georges. Principles of geostatistics. Economic geology, 58(8):1246–1266, 1963.

[3]

D. G. Krige. A statistical approach to some basic mine valuation problems on the witwatersrand. OR, 4(1):18–18, 1953. URL: http://www.jstor.org/stable/3006914 (visited on 2025-02-20).

[4]

William J. Welch, Robert J. Buck, Jerome Sacks, Henry P. Wynn, and Toby J. Mitchell. Screening, predicting, and computer experiments. Technometrics, 34(1):15–25, 1992. URL: https://www.tandfonline.com/doi/abs/10.1080/00401706.1992.10485229, arXiv:https://www.tandfonline.com/doi/pdf/10.1080/00401706.1992.10485229, doi:10.1080/00401706.1992.10485229.

[5]

Carl Edward Rasmussen and Christopher K. I. Williams. Gaussian Processes for Machine Learning. MIT Press, Cambridge, MA, 2006. ISBN 9780262182539. URL: http://www.gaussianprocess.org/gpml/.

[6]

Gregoire Allaire and Sidi Mahmoud Kaber. Numerical Linear Algebra. Texts in applied mathematics. Springer, New York, NY, January 2008.

[7]

Weiyu Liu and Stephen Batill. Gradient-enhanced response surface approximations using kriging models. In 9th AIAA/ISSMO Symposium on Multidisciplinary Analysis and Optimization. Reston, Virigina, September 2002. American Institute of Aeronautics and Astronautics.

[8]

Wataru Yamazaki, Markus Rumpfkeil, and Dimitri Mavriplis. Design optimization utilizing gradient/hessian enhanced surrogate model. In 28th AIAA Applied Aerodynamics Conference. Reston, Virigina, June 2010. American Institute of Aeronautics and Astronautics.

[9]

Selvakumar Ulaganathan, Ivo Couckuyt, Tom Dhaene, Joris Degroote, and Eric Laermans. Performance study of gradient-enhanced kriging. Eng. Comput., 32(1):15–34, January 2016.

[10]

Alexander I.J. Forrester and Andy J. Keane. Recent advances in surrogate-based optimization. Progress in Aerospace Sciences, 45(1):50–79, 2009. URL: https://www.sciencedirect.com/science/article/pii/S0376042108000766, doi:https://doi.org/10.1016/j.paerosci.2008.11.001.

[11]

Youwei He, Kuan Tan, Chunming Fu, and Jinliang Luo. An efficient gradient-enhanced kriging modeling method assisted by fast kriging for high-dimension problems. International journal of numerical methods for heat & fluid flow, 33(12):3967–3993, 2023.

[12]

Selvakumar Ulaganathan, Ivo Couckuyt, Tom Dhaene, Eric Laermans, and Joris Degroote. On the use of gradients in kriging surrogate models. In Proceedings of the Winter Simulation Conference 2014. IEEE, December 2014.

[13]

Liming Chen, Haobo Qiu, Liang Gao, Chen Jiang, and Zan Yang. A screening-based gradient-enhanced kriging modeling method for high-dimensional problems. Applied Mathematical Modelling, 69:15–31, 2019. URL: https://www.sciencedirect.com/science/article/pii/S0307904X18305900, doi:https://doi.org/10.1016/j.apm.2018.11.048.

[14]

Zhong-Hua Han, Yu Zhang, Chen-Xing Song, and Ke-Shi Zhang. Weighted gradient-enhanced kriging for high-dimensional surrogate modeling and design optimization. AIAA Journal, 55(12):4330–4346, 2017. URL: https://doi.org/10.2514/1.J055842, arXiv:https://doi.org/10.2514/1.J055842, doi:10.2514/1.J055842.

[15]

Yiming Yao, Fei Liu, and Qingfu Zhang. High-Throughput Multi-Objective bayesian optimization using gradients. In 2024 IEEE Congress on Evolutionary Computation (CEC), volume 2, 1–8. IEEE, June 2024.

[16]

Misha Padidar, Xinran Zhu, Leo Huang, Jacob R Gardner, and David Bindel. Scaling gaussian processes with derivative information using variational inference. Advances in Neural Information Processing Systems, 34:6442–6453, 2021. arXiv:2107.04061.

[17]

Haitao Liu, Jianfei Cai, and Yew-Soon Ong. Remarks on multi-output gaussian process regression. Knowledge-Based Systems, 144:102–121, 2018. URL: https://www.sciencedirect.com/science/article/pii/S0950705117306123, doi:https://doi.org/10.1016/j.knosys.2017.12.034.

[18]

B. Rakitsch, Christoph Lippert, K. Borgwardt, and Oliver Stegle. It is all in the noise: efficient multi-task gaussian process inference with structured residuals. Advances in Neural Information Processing Systems, pages, 01 2013.

[19]

Edwin V Bonilla, Kian Chai, and Christopher Williams. Multi-task gaussian process prediction. In J. Platt, D. Koller, Y. Singer, and S. Roweis, editors, Advances in Neural Information Processing Systems, volume 20. Curran Associates, Inc., 2007. URL: https://proceedings.neurips.cc/paper_files/paper/2007/file/66368270ffd51418ec58bd793f2d9b1b-Paper.pdf.

[20]

Chen Zhou Xu, Zhong Hua Han, Ke Shi Zhang, and Wen Ping Song. Improved weighted gradient-enhanced kriging model for high-dimensional aerodynamic modeling problems. In 32nd Congress of the International Council of the Aeronautical Sciences, ICAS 2021, 32nd Congress of the International Council of the Aeronautical Sciences, ICAS 2021. International Council of the Aeronautical Sciences, 2021.