Boundary Variance Inflation Causes Acquisition Bias in Gaussian Processes

# Boundary Variance Inflation Causes Acquisition Bias in Gaussian Processes
Source: [https://arxiv.org/html/2606.07561](https://arxiv.org/html/2606.07561)
,Sanna JarlRISE Research Institutes of SwedenStockholmSwedenandJens SjölundUppsala UniversityUppsalaSweden

###### Abstract\.

Gaussian processes with stationary kernels on bounded domains exhibit inflated posterior variance near the boundary\. Despite being a long\-recognized artifact in geostatistics and a source of over\-exploration in Bayesian optimization, the causes and effects of boundary\-induced acquisition bias are underexplored\. We trace the root cause to a simple geometric mechanism: the truncation of the kernel correlation neighborhood at the domain boundary creates an observation\-independent distortion that worsens with dimensionality\. We show how this distortion manifests across three acquisition classes: variance maximization concentrates selections at the corners, whereas negative integrated posterior variance and expected predictive information gain move selections inward to axis\-aligned interior shells\. These patterns arise without reference to any objective function, meaning that acquisition behavior can be dominated by kernel geometry rather than the desired task\-specific uncertainty\. To quantify this, we introduce a function\-free selection\-profile diagnostic for arbitrary acquisitions, kernels, and bounded\-domain geometries\.

Gaussian processes, posterior variance, boundary bias, Bayesian optimization, experimental design, acquisition functions, active learning

## 1\.Introduction

Gaussian processes \(GPs\) are widely used as surrogate models in Bayesian optimization \(BO\)\(Shahriariet al\.,[2016](https://arxiv.org/html/2606.07561#bib.bib10); Garnett,[2023](https://arxiv.org/html/2606.07561#bib.bib8)\)and Bayesian experimental design \(BED\)\(MacKay,[1992](https://arxiv.org/html/2606.07561#bib.bib3); Chaloner and Verdinelli,[1995](https://arxiv.org/html/2606.07561#bib.bib4); Rainforthet al\.,[2024](https://arxiv.org/html/2606.07561#bib.bib6)\)\. In both settings, the next observation is selected by optimizing acquisition functions that rely, directly or indirectly, on the GP posterior variance\. On bounded domains such as\[0,1\]D\[0,1\]^\{D\}, this variance is susceptible to a geometric artifact where the posterior variance is artificially inflated at the boundaries, regardless of the function values\.

Such boundary effects are documented in the kriging literature\(Cressie,[1993](https://arxiv.org/html/2606.07561#bib.bib21)\)\. In GP\-based sensor placement,Krauseet al\.\([2008](https://arxiv.org/html/2606.07561#bib.bib9)\)noted that entropy\-based criteria tend to cluster sensors at the boundary, “wasting” sensed information, and proposed mutual information as a remedy\.Siivolaet al\.\([2018](https://arxiv.org/html/2606.07561#bib.bib17)\)identified analogous boundary over\-exploration in BO and proposed virtual derivative sign observations as a per\-acquisition correction\. Boundary\-aware GP constructions\(Solin and Kok,[2019](https://arxiv.org/html/2606.07561#bib.bib20); Gulianet al\.,[2022](https://arxiv.org/html/2606.07561#bib.bib19)\)encode physical boundary conditions when known, while normalized convolution\(Knutsson and Westin,[1993](https://arxiv.org/html/2606.07561#bib.bib22)\)addresses the analogous loss of support in image filtering by reweighting and renormalizing convolution with a certainty field\.

Variance\-based design criteria, such as IMSPE\(Gramacy,[2020](https://arxiv.org/html/2606.07561#bib.bib14)\), mitigate certain symptoms through global integration, but the connection to the boundary mechanism is rarely made explicit\. Despite this prior work, boundary\-induced acquisition bias is seldom isolated in the development or benchmarking of acquisition functions\. The geometric mechanism is rarely separated from objective\-function behavior, the differential effect across acquisition families is not characterized, and the spatial selection patterns are not routinely measured\.

Our contribution is to isolate this variance\-inflation mechanism and characterize how a single geometric effect propagates across three variance\-driven acquisition classes\. In short, the bias stems from the kernel correlation neighborhood being truncated at the domain boundary, causing boundary points to have less in\-domain correlated volume than comparable interior points\. The effect is observation\-independent: for fixed hyperparameters, it depends on the kernel, bounded domain, and the training locations, but*not*on the observed function values\. One\-step variance maximization \(VM\), which includes the Gaussian one\-point expected information gain \(EIG\)\(Lindley,[1956](https://arxiv.org/html/2606.07561#bib.bib2)\)and the upper confidence bound \(UCB\) exploration term\(Srinivaset al\.,[2010](https://arxiv.org/html/2606.07561#bib.bib16)\), concentrates selections at corners\. Negative integrated posterior variance \(NIPV\)\(Seoet al\.,[2000](https://arxiv.org/html/2606.07561#bib.bib18); Cohnet al\.,[1996](https://arxiv.org/html/2606.07561#bib.bib7)\)and expected predictive information gain \(EPIG\)\(Bickford Smithet al\.,[2023](https://arxiv.org/html/2606.07561#bib.bib5)\)instead select axis\-aligned interior shells, with NIPV closer to the boundary than EPIG\. These patterns arise without reference to any particular objective function, showing that acquisition behavior can be dominated by kernel geometry rather than task\-specific uncertainty\.

We quantify this effect with a function\-free*selection\-profile diagnostic*, applicable to arbitrary acquisitions, kernels, and domain geometries, and validate it against sequential acquisition on representative test functions\.

## 2\.Background

We consider an unknown functionf:\[0,1\]D→ℝf\\colon\[0,1\]^\{D\}\\to\\mathbb\{R\}modeled by a GP\(Rasmussen and Williams,[2006](https://arxiv.org/html/2606.07561#bib.bib1)\),z∼𝒢​𝒫​\(m​\(𝒙\),kθ​\(𝒙,𝒙′\)\)z\\sim\\mathcal\{GP\}\(m\(\{\\bm\{x\}\}\),k\_\{\\theta\}\(\{\\bm\{x\}\},\{\\bm\{x\}\}^\{\\prime\}\)\), with noisy observationsyi=z​\(𝒙i\)\+εiy\_\{i\}=z\(\{\\bm\{x\}\}\_\{i\}\)\+\\varepsilon\_\{i\},εi∼𝒩​\(0,σn2\)\\varepsilon\_\{i\}\\sim\\mathcal\{N\}\(0,\\sigma\_\{n\}^\{2\}\)\. The posterior variance at a test point𝒙∗\{\\bm\{x\}\}\_\{\*\}and the posterior covariance between any two points are

\(1\)s2​\(𝒙∗\)=kθ​\(𝒙∗,𝒙∗\)−𝒌∗⊤​\(K\+σn2​I\)−1​𝒌∗,c​\(𝒙,𝒙′\)=kθ​\(𝒙,𝒙′\)−𝒌𝒙⊤​\(K\+σn2​I\)−1​𝒌𝒙′,s^\{2\}\(\{\\bm\{x\}\}\_\{\*\}\)=k\_\{\\theta\}\(\{\\bm\{x\}\}\_\{\*\},\{\\bm\{x\}\}\_\{\*\}\)\-\{\\bm\{k\}\}\_\{\*\}^\{\\top\}\(K\+\\sigma\_\{n\}^\{2\}I\)^\{\-1\}\{\\bm\{k\}\}\_\{\*\},\\qquad c\(\{\\bm\{x\}\},\{\\bm\{x\}\}^\{\\prime\}\)=k\_\{\\theta\}\(\{\\bm\{x\}\},\{\\bm\{x\}\}^\{\\prime\}\)\-\{\\bm\{k\}\}\_\{\\bm\{x\}\}^\{\\top\}\(K\+\\sigma\_\{n\}^\{2\}I\)^\{\-1\}\{\\bm\{k\}\}\_\{\{\\bm\{x\}\}^\{\\prime\}\},whereKi​j=kθ​\(𝒙i,𝒙j\)K\_\{ij\}=k\_\{\\theta\}\(\{\\bm\{x\}\}\_\{i\},\{\\bm\{x\}\}\_\{j\}\),𝒌𝒙=\(kθ​\(𝒙1,𝒙\),…,kθ​\(𝒙N,𝒙\)\)⊤\{\\bm\{k\}\}\_\{\\bm\{x\}\}=\(k\_\{\\theta\}\(\{\\bm\{x\}\}\_\{1\},\{\\bm\{x\}\}\),\\ldots,k\_\{\\theta\}\(\{\\bm\{x\}\}\_\{N\},\{\\bm\{x\}\}\)\)^\{\\top\}, and the variance is the diagonal cases2​\(𝒙\)=c​\(𝒙,𝒙\)s^\{2\}\(\{\\bm\{x\}\}\)=c\(\{\\bm\{x\}\},\{\\bm\{x\}\}\)\.

The class of acquisition functions whose spatial preferences are determined bys2s^\{2\}andcc\(directly or through posterior correlations\) is large, encompassing several of the most common choices\. For target points\{𝒙j∗\}j=1M\\\{\{\\bm\{x\}\}\_\{j\}^\{\*\}\\\}\_\{j=1\}^\{M\}, letcj=c​\(𝒙,𝒙j∗\)c\_\{j\}=c\(\{\\bm\{x\}\},\{\\bm\{x\}\}\_\{j\}^\{\*\}\)denote the posterior covariance between a candidate𝒙\{\\bm\{x\}\}and target𝒙j∗\{\\bm\{x\}\}\_\{j\}^\{\*\}\. The corresponding squared posterior correlation is

ρj2=cj2\(s2​\(𝒙\)\+σn2\)​\(s2​\(𝒙j∗\)\+σn2\)\.\\rho\_\{j\}^\{2\}=\\frac\{c\_\{j\}^\{2\}\}\{\(s^\{2\}\(\{\\bm\{x\}\}\)\+\\sigma\_\{n\}^\{2\}\)\(s^\{2\}\(\{\\bm\{x\}\}\_\{j\}^\{\*\}\)\+\\sigma\_\{n\}^\{2\}\)\}\.The three acquisition classes we consider, variance maximization, negative integrated posterior variance, and expected predictive information gain, are

\(2\)VM​\(𝒙\)=s2​\(𝒙\),NIPV​\(𝒙\)=∑jcj2s2​\(𝒙\)\+σn2,EPIG​\(𝒙\)=1M​∑jh​\(ρj2\),\\mathrm\{VM\}\(\{\\bm\{x\}\}\)=s^\{2\}\(\{\\bm\{x\}\}\),\\qquad\\mathrm\{NIPV\}\(\{\\bm\{x\}\}\)=\\textstyle\\sum\_\{j\}\\frac\{c\_\{j\}^\{2\}\}\{s^\{2\}\(\{\\bm\{x\}\}\)\+\\sigma\_\{n\}^\{2\}\},\\qquad\\mathrm\{EPIG\}\(\{\\bm\{x\}\}\)=\\textstyle\\frac\{1\}\{M\}\\sum\_\{j\}h\(\\rho\_\{j\}^\{2\}\),withh​\(ρ2\)=−12​log⁡\(1−ρ2\)h\(\\rho^\{2\}\)=\-\\tfrac\{1\}\{2\}\\log\(1\-\\rho^\{2\}\)\. We view variance maximization as an umbrella label that includes the GP\-UCB exploration term and GP EIG,12​log⁡\(1\+s2​\(𝒙\)/σn2\)\\tfrac\{1\}\{2\}\\log\\bigl\(1\+s^\{2\}\(\{\\bm\{x\}\}\)/\\sigma\_\{n\}^\{2\}\\bigr\), since those are monotone transforms ofs2s^\{2\}and therefore share the same argmax as pure variance maximization\(MacKay,[1992](https://arxiv.org/html/2606.07561#bib.bib3)\)\.

NIPV and EPIG both use the posterior cross\-covariancescjc\_\{j\}between candidates and targets, but aggregate them differently\. NIPV measures the total posterior variance removed by observing at𝒙\{\\bm\{x\}\}, whereas EPIG averages pairwise mutual information between candidate and target points through the convex functionhh\. We analyze EPIG in the prediction\-oriented setting where targets span the full domain, so all three methods share the same candidate space and can be compared directly\. In BO, full UCBμ​\(𝒙\)\+β​σ​\(𝒙\)\\mu\(\{\\bm\{x\}\}\)\+\\beta\\sigma\(\{\\bm\{x\}\}\)depends on𝒚\{\\bm\{y\}\}throughμ\\mu, but its exploration component does not\. When exploration dominates early in optimization, this variance bias controls where observations are placed\.

## 3\.Boundary Variance Inflation

The posterior variances2​\(𝒙∗\)s^\{2\}\(\{\\bm\{x\}\}\_\{\*\}\)is the prior variance minus a reduction determined by the covariance between𝒙∗\{\\bm\{x\}\}\_\{\*\}and the training inputs\. For stationary kernels, this covariance decays with distance at a rate determined by the lengthscaleℓ\\ell\. On a bounded domain, the correlation neighborhood of a boundary point is truncated: training data can contribute from the interior side but not from beyond the boundary\.

###### Definition 3\.0 \(Boundary distance\)\.

For𝒙∈\[0,1\]D\{\\bm\{x\}\}\\in\[0,1\]^\{D\}, the boundary distance isd∂​\(𝒙\)=mind⁡min⁡\(xd,1−xd\)d\_\{\\partial\}\(\{\\bm\{x\}\}\)=\\min\_\{d\}\\min\(x\_\{d\},1\-x\_\{d\}\), the minimum distance from𝒙\{\\bm\{x\}\}to any face of the unit hypercube\. This ranges from0on any face to1/21/2at the cube’s center\.

An interior point withd∂​\(𝒙\)≫ℓd\_\{\\partial\}\(\{\\bm\{x\}\}\)\\gg\\ellhas an approximately complete correlation neighborhood within the domain; a boundary point withd∂​\(𝒙\)≪ℓd\_\{\\partial\}\(\{\\bm\{x\}\}\)\\ll\\ellhas this neighborhood truncated by one or more faces\. On average, the boundary point has higher posterior variance than an interior point with comparable local training density\. The effect compounds with dimension: the boundary truncates the correlation neighborhood independently in each coordinate, so corner points lose correlated volume along allDDfaces simultaneously\. For typicalNNandDD, this systematically elevates posterior variance at corners and edges relative to the interior\. Because this inflated variance enters the acquisition functions through different functional expressions, the resulting artifacts differ\.

![Refer to caption](https://arxiv.org/html/2606.07561v1/x1.png)Figure 1\.Acquisition surfaces and argmax\-location densities on\[0,1\]2\[0,1\]^\{2\}under an isotropic Matérn\-5/2\\nicefrac\{\{5\}\}\{\{2\}\}kernel \(ℓ=0\.2\\ell=0\.2,N=15N=15\)\.Left:Acquisition value, individually normalized, for one fixed Sobol training design;★\\bigstarmarks the argmax\.Right:D4\-symmetrized argmax\-location density over10 00010\\,000random training sets\. VM concentrates at the corners, NIPV forms a near\-boundary square frame, and EPIG forms a frame further inside\. These are 2D instances of the axis\-aligned interior shells discussed throughout, determined by kernel geometry and training locations rather than observed function values\.Two\-by\-three grid of acquisition surfaces \(left\) and argmax densities \(right\) on the unit square for VM, NIPV, and EPIG\. VM peaks at corners, NIPV along edges away from corners, EPIG in the interior\.VM inherits the artifact directly, concentrating observations at the boundary because the acquisition is monotone ins2​\(𝒙\)s^\{2\}\(\{\\bm\{x\}\}\)\. NIPV partially corrects it: the high local variances2​\(𝒙\)s^\{2\}\(\{\\bm\{x\}\}\)at a boundary candidate inflates the denominator, while the smaller cross\-covariancescjc\_\{j\}with the target set depress the numerator\. The maximum sits on an axis\-aligned interior shell just inside the boundary\. EPIG shifts selection further inward through the soft\-threshold structure ofhh, which is steep nearρ2=1\\rho^\{2\}=1and shallow near zero: each target contributes substantially only if the candidate has high posterior correlation with it\. An interior candidate achieves this with targets in all directions; a boundary candidate’s truncated correlations mean that most distant targets contribute little to the average\. Consequently, EPIG’s maximum sits on a shell at a largerd∂d\_\{\\partial\}than NIPV’s\.

Under fixed hyperparameters, none of these spatial preferences responds to the observed function values \([Remark2\.1](https://arxiv.org/html/2606.07561#S2.Thmremark1)\); in particular, none can, by itself, identify regions where the realized function varies most\. The boundary bias is therefore a property of the variance field itself; the acquisition function determines only how that bias is translated into a placement preference\.

[Figure1](https://arxiv.org/html/2606.07561#S3.F1)illustrates the effect in two dimensions\. The left panels show one representative acquisition surface for each method, computed for a fixed Sobol design under an isotropic Matérn\-5/2\\nicefrac\{\{5\}\}\{\{2\}\}kernel \(ℓ=0\.2\\ell=0\.2,N=15N=15\)\. VM peaks sharply at a domain corner \(★\\bigstar\); NIPV peaks near the boundary but away from the corner; EPIG peaks in the interior\. The right panels show the argmax\-location density over 10 000 random training designs\. Together, these show the bias reflects kernel geometry and training locations, not observed function values\.

## 4\.The Selection Profile Diagnostic

The diagnostic tracks where the argmax occurs, rather than the mean acquisition value at each boundary distance\. This matters because sequential design selects maximizers, and mean profiles can be nearly flat even when the argmax repeatedly falls in specific boundary\-distance shells\. We define the*selection profile*as the empirical distribution ofd∂​\(𝒙next\)d\_\{\\partial\}\(\{\\bm\{x\}\}\_\{\\mathrm\{next\}\}\)over many independent training sets\. This gives a one\-dimensional summary in any dimension\.

We run a*controlled diagnostic*: for fixedDD,NN, and kernel hyperparameters, we drawS=1000S=1000random training sets on\[0,1\]D\[0,1\]^\{D\}, compute the GP posterior under a Matérn\-5/2\\nicefrac\{\{5\}\}\{\{2\}\}kernel, and record the boundary\-distance bin of the acquisition argmax\. To separate acquisition preference from cube geometry, candidates are stratified across boundary\-distance shells and then reweighted to recover the uniform\-sampling distribution\.

![Refer to caption](https://arxiv.org/html/2606.07561v1/x2.png)Figure 2\.Selection profiles under a Matérn\-5/2\\nicefrac\{\{5\}\}\{\{2\}\}kernel withN=50N=50training points \(1000 seeds; shaded bands are 95% confidence intervals\)\. Each curve shows the selection probability of each boundary\-distance bin after reweighting by shell volume; the dotted baseline is the geometric no\-preference distribution\. Lengthscalesℓ=0\.20,0\.30,0\.40,0\.55\\ell=0\.20,0\.30,0\.40,0\.55forD=2,3,4,6D=2,3,4,6are chosen as the smallest sweep value above the nearest\-neighbor distanceN−1/DN^\{\-1/D\};[SectionC\.2](https://arxiv.org/html/2606.07561#A3.SS2)shows the bias is qualitatively robust across this range\.Four\-panel line plot of selection probability versus boundary distance, one panel per dimension $D\\in\\\{2,3,4,6\\\}$\. VM concentrates sharply at the boundary in every panel; NIPV and EPIG peak at interior shells, with NIPV’s closer to the boundary and the two converging as $D$ grows\.[Figure2](https://arxiv.org/html/2606.07561#S4.F2)shows the result\. The 2D pattern of[Figure1](https://arxiv.org/html/2606.07561#S3.F1)persists across dimensions: VM stays at the boundary, while NIPV and EPIG settle onto interior shells whose locations depend onℓ\\ellandDD\. In lowDD, NIPV lies closer to the boundary; the shells converge asDDgrows\. Mean acquisition\-value profiles miss these patterns because the argmax is governed by the upper tail, not the mean \([SectionB\.1](https://arxiv.org/html/2606.07561#A2.SS1)\)\.

The controlled diagnostic also matches sequential acquisition on representative test functions \([SectionC\.1](https://arxiv.org/html/2606.07561#A3.SS1)\), confirming that the bias reflects kernel geometry rather than test\-function properties\.

## 5\.Discussion

Boundary variance inflation affects both BO and BED because it enters through the GP variance field\. In BO, this is not specific to UCB: any acquisition whose exploratory behavior is driven by posterior uncertainty can inherit the bias when exploration dominates\. We illustrate this with UCB, whose exploration term is explicit, and give a one\-step BO example in[AppendixD](https://arxiv.org/html/2606.07561#A4)\. Related effects can also arise in safe BO\(Suiet al\.,[2015](https://arxiv.org/html/2606.07561#bib.bib23)\): inflated boundary uncertainty widens confidence bands and can make boundary regions harder to certify as safe\. In information\-based design, the same effect can favor geometrically determined locations over the most informative ones\. Thus, the issue is not a peculiarity of one acquisition function, but a shared consequence of using a geometrically distorted variance field\.

A Neumann \(mirror\-image\) kernel mitigates the root cause by completing the truncated correlation neighborhood with reflected images \([AppendixE](https://arxiv.org/html/2606.07561#A5)\)\. As shown in[Figure8](https://arxiv.org/html/2606.07561#A5.F8), it reduces VM’s boundary concentration, removing7%7\\%of the excess over the geometric baseline atD=2D=2and28%28\\%atD=6D=6\. The correction is partial: the image method changes the covariance geometry without adding independent observations beyond the domain\. More generally, reducing geometric nonuniformity in the variance field would decrease the shared source of bias across acquisition functions\.

Our analysis isolates the variance\-field mechanism under fixed hyperparameters; in full sequential workflows, hyperparameter learning and posterior\-mean terms can attenuate or amplify it\. The selection\-profile diagnostic provides a function\-free way to check arbitrary acquisitions, kernels, and domain geometries for unintended spatial preferences\. Making these preferences measurable can help guide the design of acquisition functions and kernels with more reliable sampling behavior\.

## Code Availability

## Author Contributions

Maria Bånkestad: Conceptualization; Methodology; Software; Formal Analysis; Investigation; Visualization; Writing\-Original Draft; Writing\-Review & Editing\.Sanna Jarl: Investigation; Writing\-Review & Editing\.Jens Sjölund: Conceptualization; Methodology; Writing\-Review & Editing\.

## References

- F\. Bickford Smith, A\. Kirsch, S\. Farquhar, Y\. Gal, A\. Foster, and T\. Rainforth \(2023\)Prediction\-oriented Bayesian active learning\.InProceedings of The 26th International Conference on Artificial Intelligence and Statistics,Proceedings of Machine Learning Research, Vol\.206,Valencia, Spain,pp\. 7331–7348\.Cited by:[§1](https://arxiv.org/html/2606.07561#S1.p4.1)\.
- K\. Chaloner and I\. Verdinelli \(1995\)Bayesian experimental design: a review\.Statistical Science10\(3\),pp\. 273–304\.Cited by:[§1](https://arxiv.org/html/2606.07561#S1.p1.1)\.
- D\. A\. Cohn, Z\. Ghahramani, and M\. I\. Jordan \(1996\)Active learning with statistical models\.Journal of Artificial Intelligence Research4,pp\. 129–145\.Cited by:[§1](https://arxiv.org/html/2606.07561#S1.p4.1)\.
- N\. A\. C\. Cressie \(1993\)Statistics for spatial data\.Revised edition,Wiley,New York\.Cited by:[§1](https://arxiv.org/html/2606.07561#S1.p2.1)\.
- R\. Garnett \(2023\)Bayesian optimization\.Cambridge University Press,Cambridge, UK\.Cited by:[§1](https://arxiv.org/html/2606.07561#S1.p1.1)\.
- R\. B\. Gramacy \(2020\)Surrogates: Gaussian process modeling, design, and optimization for the applied sciences\.Chapman and Hall/CRC,Boca Raton, FL, USA\.Cited by:[§1](https://arxiv.org/html/2606.07561#S1.p3.1)\.
- M\. Gulian, A\. Frankel, and L\. Swiler \(2022\)Gaussian process regression constrained by boundary value problems\.Computer Methods in Applied Mechanics and Engineering388,pp\. 114117\.External Links:[Document](https://dx.doi.org/10.1016/j.cma.2021.114117)Cited by:[§1](https://arxiv.org/html/2606.07561#S1.p2.1)\.
- C\. Hvarfner, E\. O\. Hellsten, and L\. Nardi \(2024\)Vanilla Bayesian optimization performs great in high dimensions\.InProceedings of the 41st International Conference on Machine Learning,Proceedings of Machine Learning Research, Vol\.235,Vienna, Austria,pp\. 20793–20817\.Cited by:[§C\.1](https://arxiv.org/html/2606.07561#A3.SS1.p1.8)\.
- H\. Knutsson and C\. Westin \(1993\)Normalized and differential convolution: methods for interpolation and filtering of incomplete and uncertain data\.InProceedings of the IEEE Conference on Computer Vision and Pattern Recognition \(CVPR\),New York, NY, USA,pp\. 515–523\.Cited by:[§1](https://arxiv.org/html/2606.07561#S1.p2.1)\.
- A\. Krause, A\. Singh, and C\. Guestrin \(2008\)Near\-optimal sensor placements in Gaussian processes: theory, efficient algorithms and empirical studies\.Journal of Machine Learning Research9,pp\. 235–284\.Cited by:[§1](https://arxiv.org/html/2606.07561#S1.p2.1)\.
- D\. V\. Lindley \(1956\)On a measure of the information provided by an experiment\.The Annals of Mathematical Statistics27\(4\),pp\. 986–1005\.Cited by:[§1](https://arxiv.org/html/2606.07561#S1.p4.1)\.
- D\. J\. C\. MacKay \(1992\)Information\-based objective functions for active data selection\.Neural Computation4\(4\),pp\. 590–604\.External Links:[Document](https://dx.doi.org/10.1162/neco.1992.4.4.590)Cited by:[§1](https://arxiv.org/html/2606.07561#S1.p1.1),[§2](https://arxiv.org/html/2606.07561#S2.p2.9)\.
- T\. Rainforth, A\. Foster, D\. R\. Ivanova, and F\. Bickford Smith \(2024\)Modern Bayesian experimental design\.Statistical Science39\(1\),pp\. 100–114\.Cited by:[§1](https://arxiv.org/html/2606.07561#S1.p1.1)\.
- C\. E\. Rasmussen and C\. K\. I\. Williams \(2006\)Gaussian processes for machine learning\.MIT Press,Cambridge, MA, USA\.Cited by:[§2](https://arxiv.org/html/2606.07561#S2.p1.5)\.
- S\. Seo, M\. Wallat, T\. Graepel, and K\. Obermayer \(2000\)Gaussian process regression: active data selection and test point rejection\.InProceedings of the IEEE\-INNS\-ENNS International Joint Conference on Neural Networks \(IJCNN\),Vol\.3,Como, Italy,pp\. 241–246\.External Links:[Document](https://dx.doi.org/10.1109/IJCNN.2000.861310)Cited by:[§1](https://arxiv.org/html/2606.07561#S1.p4.1)\.
- B\. Shahriari, K\. Swersky, Z\. Wang, R\. P\. Adams, and N\. De Freitas \(2016\)Taking the human out of the loop: a review of Bayesian optimization\.Proceedings of the IEEE104\(1\),pp\. 148–175\.Cited by:[§1](https://arxiv.org/html/2606.07561#S1.p1.1)\.
- E\. Siivola, A\. Vehtari, J\. Vanhatalo, J\. González, and M\. R\. Andersen \(2018\)Correcting boundary over\-exploration deficiencies in Bayesian optimization with virtual derivative sign observations\.In2018 IEEE 28th International Workshop on Machine Learning for Signal Processing \(MLSP\),Aalborg, Denmark,pp\. 1–6\.External Links:[Document](https://dx.doi.org/10.1109/MLSP.2018.8516936)Cited by:[§1](https://arxiv.org/html/2606.07561#S1.p2.1)\.
- J\. Snoek, K\. Swersky, R\. Zemel, and R\. Adams \(2014\)Input warping for Bayesian optimization of non\-stationary functions\.InProceedings of the 31st International Conference on Machine Learning,Proceedings of Machine Learning Research, Vol\.32,Beijing, China,pp\. 1674–1682\.Cited by:[Appendix E](https://arxiv.org/html/2606.07561#A5.p7.1)\.
- A\. Solin and M\. Kok \(2019\)Know your boundaries: constraining Gaussian processes by variational harmonic features\.InProceedings of the Twenty\-Second International Conference on Artificial Intelligence and Statistics,Proceedings of Machine Learning Research, Vol\.89,Naha, Okinawa, Japan,pp\. 2193–2202\.Cited by:[§1](https://arxiv.org/html/2606.07561#S1.p2.1)\.
- N\. Srinivas, A\. Krause, S\. M\. Kakade, and M\. Seeger \(2010\)Gaussian process optimization in the bandit setting: no regret and experimental design\.InInternational Conference on Machine Learning,Haifa, Israel,pp\. 1015–1022\.Cited by:[Appendix D](https://arxiv.org/html/2606.07561#A4.p2.13),[§1](https://arxiv.org/html/2606.07561#S1.p4.1)\.
- Y\. Sui, A\. Gotovos, J\. Burdick, and A\. Krause \(2015\)Safe exploration for optimization with Gaussian processes\.InProceedings of the 32nd International Conference on Machine Learning,Proceedings of Machine Learning Research, Vol\.37,Lille, France,pp\. 997–1005\.Cited by:[§5](https://arxiv.org/html/2606.07561#S5.p1.1)\.
- S\. Surjanovic and D\. Bingham \(2013\)Virtual library of simulation experiments: test functions and datasets\.Note:Accessed: 2026\-05\-25External Links:[Link](https://www.sfu.ca/~ssurjano/optimization.html)Cited by:[§C\.1](https://arxiv.org/html/2606.07561#A3.SS1.p1.8),[Appendix D](https://arxiv.org/html/2606.07561#A4.p2.13)\.
- R\. Vershynin \(2018\)High\-dimensional probability: an introduction with applications in data science\.Cambridge University Press,Cambridge, UK\.Cited by:[§B\.3](https://arxiv.org/html/2606.07561#A2.SS3.p1.5)\.

## Appendix AAdditional Theory and Mechanism

### A\.1\.Gaussian Process Posterior

Given observations𝒟=\{X,𝒚\}\\mathcal\{D\}=\\\{X,\{\\bm\{y\}\}\\\}, the GP predictive posterior at a test point𝒙∗\{\\bm\{x\}\}\_\{\*\}is Gaussian with mean and variance

\(3\)μ​\(𝒙∗\)\\displaystyle\\mu\(\{\\bm\{x\}\}\_\{\*\}\)=m​\(𝒙∗\)\+𝒌∗⊤​\(K\+σn2​I\)−1​\(𝒚−m​\(X\)\),\\displaystyle=m\(\{\\bm\{x\}\}\_\{\*\}\)\+\{\\bm\{k\}\}\_\{\*\}^\{\\top\}\(K\+\\sigma\_\{n\}^\{2\}I\)^\{\-1\}\(\{\\bm\{y\}\}\-m\(X\)\),\(4\)s2​\(𝒙∗\)\\displaystyle s^\{2\}\(\{\\bm\{x\}\}\_\{\*\}\)=kθ​\(𝒙∗,𝒙∗\)−𝒌∗⊤​\(K\+σn2​I\)−1​𝒌∗,\\displaystyle=k\_\{\\theta\}\(\{\\bm\{x\}\}\_\{\*\},\{\\bm\{x\}\}\_\{\*\}\)\-\{\\bm\{k\}\}\_\{\*\}^\{\\top\}\(K\+\\sigma\_\{n\}^\{2\}I\)^\{\-1\}\{\\bm\{k\}\}\_\{\*\},whereKi​j=kθ​\(𝒙i,𝒙j\)K\_\{ij\}=k\_\{\\theta\}\(\{\\bm\{x\}\}\_\{i\},\{\\bm\{x\}\}\_\{j\}\)and𝒌𝒙=\(kθ​\(𝒙1,𝒙\),…,kθ​\(𝒙N,𝒙\)\)⊤\{\\bm\{k\}\}\_\{\\bm\{x\}\}=\(k\_\{\\theta\}\(\{\\bm\{x\}\}\_\{1\},\{\\bm\{x\}\}\),\\ldots,k\_\{\\theta\}\(\{\\bm\{x\}\}\_\{N\},\{\\bm\{x\}\}\)\)^\{\\top\}is the vector of kernel evaluations between the training points and𝒙\{\\bm\{x\}\}\. The posterior covariance between two test points is

\(5\)c​\(𝒙,𝒙′\)=kθ​\(𝒙,𝒙′\)−𝒌𝒙⊤​\(K\+σn2​I\)−1​𝒌𝒙′,c\(\{\\bm\{x\}\},\{\\bm\{x\}\}^\{\\prime\}\)=k\_\{\\theta\}\(\{\\bm\{x\}\},\{\\bm\{x\}\}^\{\\prime\}\)\-\{\\bm\{k\}\}\_\{\\bm\{x\}\}^\{\\top\}\(K\+\\sigma\_\{n\}^\{2\}I\)^\{\-1\}\{\\bm\{k\}\}\_\{\{\\bm\{x\}\}^\{\\prime\}\},withs2​\(𝒙\)=c​\(𝒙,𝒙\)s^\{2\}\(\{\\bm\{x\}\}\)=c\(\{\\bm\{x\}\},\{\\bm\{x\}\}\)as the diagonal case\. The posterior meanμ​\(𝒙\)\\mu\(\{\\bm\{x\}\}\)depends on the observed values𝒚\{\\bm\{y\}\}, but neithers2s^\{2\}norccdoes\. This is immediate from the expressions above: both are determined entirely byXX,kθk\_\{\\theta\}, andσn2\\sigma\_\{n\}^\{2\}\.

### A\.2\.One\-Dimensional Variance Reduction Example

To make the truncated\-neighborhood argument concrete, consider a Matérn\-5/2\\nicefrac\{\{5\}\}\{\{2\}\}kernel in one dimension\. For a single training point𝒙1\{\\bm\{x\}\}\_\{1\}at distancehhfrom𝒙∗\{\\bm\{x\}\}\_\{\*\}, the variance reduction is

Δ​s2=kθ​\(𝒙∗,𝒙1\)2kθ​\(𝒙1,𝒙1\)\+σn2,\\Delta s^\{2\}=\\frac\{k\_\{\\theta\}\(\{\\bm\{x\}\}\_\{\*\},\{\\bm\{x\}\}\_\{1\}\)^\{2\}\}\{k\_\{\\theta\}\(\{\\bm\{x\}\}\_\{1\},\{\\bm\{x\}\}\_\{1\}\)\+\\sigma\_\{n\}^\{2\}\},which depends onh/ℓh/\\ellbut not on the position of𝒙∗\{\\bm\{x\}\}\_\{\*\}relative to the boundary\. With multiple training points, the reduction is governed by the full covariance vector and kernel matrix, but the same intuition holds: fewer in\-domain correlated observations are available near the boundary\. A boundary point, therefore, receives less variance reduction on average, producing the inflation visible in[Figure3](https://arxiv.org/html/2606.07561#A1.F3)\.

![Refer to caption](https://arxiv.org/html/2606.07561v1/x3.png)Figure 3\.Posterior variances2​\(x\)s^\{2\}\(x\)for a Matérn\-5/2\\nicefrac\{\{5\}\}\{\{2\}\}kernel \(ℓ=0\.12\\ell=0\.12,σn2=0\.003\\sigma\_\{n\}^\{2\}=0\.003\) given five evenly spaced training points on\[0,1\]\[0,1\]\. Interior variance dips toward zero at each training point and rises moderately between them\. Near the boundary \(x=0x=0andx=1x=1, shaded\), the kernel’s correlation neighborhood is truncated and the variance climbs back toward the prior level \(dashed\): the geometric artifact that drives the bias\.One\-dimensional plot of GP posterior variance versus position on the unit interval\. Variance dips toward zero at five interior training points and rises sharply toward the prior level near the two boundaries\.
### A\.3\.Detailed Inheritance Analysis

#### Pointwise variance criteria: boundary concentration\.

For a GP with Gaussian likelihood,EIG​\(𝒙\)=12​log⁡\(1\+s2​\(𝒙\)/σn2\)\\mathrm\{EIG\}\(\{\\bm\{x\}\}\)=\\tfrac\{1\}\{2\}\\log\\bigl\(1\+s^\{2\}\(\{\\bm\{x\}\}\)/\\sigma\_\{n\}^\{2\}\\bigr\)is a monotone transform ofs2​\(𝒙\)s^\{2\}\(\{\\bm\{x\}\}\)\. The boundary variance inflation therefore passes through directly: under the conditions of[Section3](https://arxiv.org/html/2606.07561#S3),𝔼X​\[EIG​\(𝒙bnd\)\]\>𝔼X​\[EIG​\(𝒙int\)\]\\mathbb\{E\}\_\{X\}\[\\mathrm\{EIG\}\(\{\\bm\{x\}\}\_\{\\mathrm\{bnd\}\}\)\]\>\\mathbb\{E\}\_\{X\}\[\\mathrm\{EIG\}\(\{\\bm\{x\}\}\_\{\\mathrm\{int\}\}\)\], so the EIG argmax tends toward boundary locations\. The same holds for UCB’s exploration termβ​σ​\(𝒙\)\\beta\\sigma\(\{\\bm\{x\}\}\)and for pure variance maximization, since all three share the same variance\-driven argmax\. Pointwise variance criteria therefore allocate observations according to kernel geometry rather than task\-specific uncertainty\.

#### NIPV: partial correction through global integration\.

NIPV evaluates each candidate by the extent to which it reduces variance across the entire domain\. A boundary candidate𝒙\{\\bm\{x\}\}has high local variances2​\(𝒙\)s^\{2\}\(\{\\bm\{x\}\}\), but its posterior covariancesc​\(𝒙,𝒙j∗\)c\(\{\\bm\{x\}\},\{\\bm\{x\}\}\_\{j\}^\{\*\}\)with distant target points are truncated by the boundary, limiting the numerator∑jcj2\\sum\_\{j\}c\_\{j\}^\{2\}\. An interior candidate has lower local variance but correlates symmetrically with a larger volume of the domain\. This trade\-off displaces NIPV’s argmax onto an interior shell just inside the boundary\. The shell’s location depends on the lengthscale and dimension\.

#### EPIG: a characteristic interior shell\.

EPIG averages pairwise mutual information termsh​\(ρj2\)=−12​log⁡\(1−ρj2\)h\(\\rho\_\{j\}^\{2\}\)=\-\\frac\{1\}\{2\}\\log\(1\-\\rho\_\{j\}^\{2\}\)over target points\. The functionhhis convex and steep nearρ2=1\\rho^\{2\}=1\. This nonlinearity creates a qualitative difference from both EIG and NIPV:

1. \(i\)For a boundary candidate withd∂​\(𝒙\)≪ℓd\_\{\\partial\}\(\{\\bm\{x\}\}\)\\ll\\ell, nearby targets on the interior side have highρj2\\rho\_\{j\}^\{2\}and contribute large EPIG terms\. But targets on the far side of the domain have lower posterior correlation with𝒙\{\\bm\{x\}\}than they would with an interior candidate, because the boundary truncates the propagation of correlations\. The average over all targets is pulled down by these many poorly correlated distant targets\.
2. \(ii\)For an interior candidate withd∂​\(𝒙\)≫ℓd\_\{\\partial\}\(\{\\bm\{x\}\}\)\\gg\\ell, the correlation field is symmetric and extends in all directions\. Targets at moderate distance contribute a substantialρj2\\rho\_\{j\}^\{2\}from every direction, and the average is not truncated\.

The averaging acts as a spatial regularizer, penalizing boundary candidates\. But it does not produce spatially unbiased sampling\. The EPIG argmax concentrates at a characteristic interior shell at a boundary distance that depends onℓ\\ellandDD\. This shell arises because EPIG favors candidates with symmetric, far\-reaching correlation neighborhoods, and there is a specific distance from the boundary where this symmetry is maximized relative to the competing effect of elevated prior variance\.

#### NIPV\-EPIG convergence asDDgrows\.

NIPV and EPIG share the same cross\-covariance informationcjc\_\{j\}but aggregate it differently: NIPV linearly, EPIG through the saturatingh​\(ρj2\)h\(\\rho\_\{j\}^\{2\}\)\. Forρj2≪1\\rho\_\{j\}^\{2\}\\ll 1, the Taylor expansion givesh​\(ρj2\)≈ρj2/2h\(\\rho\_\{j\}^\{2\}\)\\approx\\rho\_\{j\}^\{2\}/2, and substituting yields

EPIG​\(𝒙\)≈12​M​\(s2​\(𝒙\)\+σn2\)​∑jcj2s2​\(𝒙j∗\)\+σn2,\\mathrm\{EPIG\}\(\{\\bm\{x\}\}\)\\;\\approx\\;\\frac\{1\}\{2M\\,\(s^\{2\}\(\{\\bm\{x\}\}\)\+\\sigma\_\{n\}^\{2\}\)\}\\sum\_\{j\}\\frac\{c\_\{j\}^\{2\}\}\{s^\{2\}\(\{\\bm\{x\}\}\_\{j\}^\{\*\}\)\+\\sigma\_\{n\}^\{2\}\},which is NIPV with an additional per\-target reweighting1/\(s2​\(𝒙j∗\)\+σn2\)1/\(s^\{2\}\(\{\\bm\{x\}\}\_\{j\}^\{\*\}\)\+\\sigma\_\{n\}^\{2\}\)\. If target variances are roughly uniform across the target set, the reweighting collapses to a constant andEPIG∝NIPV\\mathrm\{EPIG\}\\propto\\mathrm\{NIPV\}– the two share an argmax\. Two effects pushρj2\\rho\_\{j\}^\{2\}into this small\-correlation regime asDDgrows: withNNfixed andℓ∼N−1/D\\ell\\sim N^\{\-1/D\}, sparse training keeps target variances close to the prior level and enlarges the EPIG denominator; and the product kernelk​\(𝒙,𝒙j∗\)=∏dκ​\(xd−xj,d∗\)k\(\{\\bm\{x\}\},\{\\bm\{x\}\}\_\{j\}^\{\*\}\)=\\prod\_\{d\}\\kappa\(x\_\{d\}\-x\_\{j,d\}^\{\*\}\)attenuates sharply as soon as any coordinate distance exceedsℓ\\ell, so fewer candidate\-target pairs reach high posterior correlation\. The convergence is gradual: profiles are noticeably similar fromD=4D=4and essentially coincide atD=6D=6\(right panel of[Figure2](https://arxiv.org/html/2606.07561#S4.F2)\)\.

## Appendix BSelection\-Profile Diagnostic Details

### B\.1\.Why the Argmax Matters, Not the Mean

A natural first attempt at diagnosing spatial bias is to plot the mean acquisition value as a function of boundary distance:𝔼​\[α​\(𝒙\)∣d∂​\(𝒙\)∈Bb\]\\mathbb\{E\}\[\\alpha\(\{\\bm\{x\}\}\)\\mid d\_\{\\partial\}\(\{\\bm\{x\}\}\)\\in B\_\{b\}\]\. This turns out to be misleading\. In our experiments, EPIG’s mean profile is nearly flat across boundary\-distance bins, suggesting a lack of spatial preference\. Yet the argmax of EPIG consistently falls in interior bins\. The reason is that sequential design selects the maximizer𝒙next=arg​max𝒙∈𝒳c⁡α​\(𝒙\)\{\\bm\{x\}\}\_\{\\mathrm\{next\}\}=\\operatorname\*\{arg\\,max\}\_\{\{\\bm\{x\}\}\\in\\mathcal\{X\}\_\{c\}\}\\alpha\(\{\\bm\{x\}\}\), so what determines sampling behavior is the upper tail of the acquisition distribution within each bin, not the conditional mean\. A bin with slightly higher variance in its acquisition values will more often produce the global maximum, even if its mean is the same as that of other bins\. This is why the selection profile records argmax locations rather than mean acquisition values\.

A related concern is that, in practice, BO uses approximate acquisition optimization rather than the exact grid argmax, which could, in principle, dilute the bias by missing the exact maximizer\. Two observations argue against this\. First, the bias concerns*where*the high\-acquisition regions sit: any optimizer that climbs the acquisition surface is pulled into the boundary\-inflated regions, regardless of whether it recovers the exact argmax\. Second, the sequential validation of[Figure4](https://arxiv.org/html/2606.07561#A3.F4)already uses MAP\-fitted GPs with realistic acquisition optimization on real test functions, and the spatial bias persists: it is not an artifact of exact\-argmax recovery on a candidate grid\.

### B\.2\.Controlled Diagnostic Setup

We fix the dimensionDD, training set sizeNN, and kernel hyperparameters\(ℓ,σn2\)\(\\ell,\\sigma\_\{n\}^\{2\}\), then repeat the followingS=1000S=1000times\. DrawNNtraining points uniformly on\[0,1\]D\[0,1\]^\{D\}and compute the GP posterior under a product Matérn\-5/2\\nicefrac\{\{5\}\}\{\{2\}\}kernel\. Evaluate each acquisition function on a candidate set𝒳c\\mathcal\{X\}\_\{c\}and record which boundary\-distance bin contains the argmax\. The selection probability for binbbis

\(6\)p^b=1S​∑s=1S𝟏​\[d∂​\(𝒙next\(s\)\)∈Bb\]\.\\hat\{p\}\_\{b\}\\;=\\;\\frac\{1\}\{S\}\\sum\_\{s=1\}^\{S\}\\mathbf\{1\}\\\!\\bigl\[d\_\{\\partial\}\(\{\\bm\{x\}\}\_\{\\mathrm\{next\}\}^\{\(s\)\}\)\\in B\_\{b\}\\bigr\]\.
For the headline figures, the lengthscale at eachDDis chosen as the smallest value in the sweep grid of[SectionC\.2](https://arxiv.org/html/2606.07561#A3.SS2)above the typical nearest\-neighbor distanceN−1/DN^\{\-1/D\}ofN=50N=50uniform training points; this places the kernel range2​ℓ2\\ellcomfortably above the inter\-point spacing, the regime in which the GP meaningfully interpolates between training points\. The specific settings areN=50N=50throughout, withℓ=0\.20\\ell=0\.20forD=2D=2,ℓ=0\.30\\ell=0\.30forD=3D=3,ℓ=0\.40\\ell=0\.40forD=4D=4, andℓ=0\.55\\ell=0\.55forD=6D=6; observation noiseσn2=0\.003\\sigma\_\{n\}^\{2\}=0\.003\. For the 2D surface figure \([Figure1](https://arxiv.org/html/2606.07561#S3.F1)\), we use an isotropic Matérn\-5/2\\nicefrac\{\{5\}\}\{\{2\}\}kernel withℓ=0\.2\\ell=0\.2,N=15N=15,σn2=0\.003\\sigma\_\{n\}^\{2\}=0\.003, evaluated on a200×200200\\times 200grid with 2048 Sobol integration points; the argmax density aggregates 10 000 random training sets on a128×128128\\times 128candidate grid\.

### B\.3\.Stratification and Volume Reweighting

The unit hypercube\[0,1\]D\[0,1\]^\{D\}has far more volume near its boundary than in its interior\(Vershynin,[2018](https://arxiv.org/html/2606.07561#bib.bib11)\)\. Under uniform sampling, the density ofd∂d\_\{\\partial\}isf​\(t\)=2​D​\(1−2​t\)D−1f\(t\)=2D\(1\-2t\)^\{D\-1\}, which is sharply concentrated neart=0t=0for moderateDD\. If candidates are drawn uniformly, the selection profile mixes two effects: the acquisition function’s preference and the cube’s geometric volume\.

To separate them, we build a stratified candidate set in which each boundary\-distance bin contains the same number of candidates\. This gives each shell equal opportunity to contain the argmax and reveals the acquisition’s intrinsic preference overd∂d\_\{\\partial\}\. For plots that represent uniform candidate availability in the cube, we then reweight the stratified bin probabilities by the shell volumes

wb∝\(1−2​ab\)D−\(1−2​bb\)D,w\_\{b\}\\propto\(1\-2a\_\{b\}\)^\{D\}\-\(1\-2b\_\{b\}\)^\{D\},whereaba\_\{b\}andbbb\_\{b\}are the bin edges\. The dashed baseline is the corresponding no\-preference distribution: the distribution obtained by an acquisition function that is equally likely to select any uniformly available candidate\.

## Appendix CRobustness and Sequential Validation

### C\.1\.Sequential Validation

To confirm that the controlled diagnostic predicts real sequential behavior, we run acquisition on three test functions: Peaks \(MATLAB\) \(D=2D=2,N0=35N\_\{0\}=35\), Smooth\-box, a product of one\-dimensional smooth step functions \(D=3D=3,N0=40N\_\{0\}=40\), and Hartmann\(Surjanovic and Bingham,[2013](https://arxiv.org/html/2606.07561#bib.bib12)\)\(D=4D=4,N0=50N\_\{0\}=50\)\. Starting fromN0N\_\{0\}initial points, we fit a GP with MAP hyperparameters under a dimension\-scaled log\-normal lengthscale prior\(Hvarfneret al\.,[2024](https://arxiv.org/html/2606.07561#bib.bib13)\)and select one point via each acquisition function\. Using 1000 independent seeds, we histogramd∂​\(𝒙next\)d\_\{\\partial\}\(\{\\bm\{x\}\}\_\{\\mathrm\{next\}\}\)and compare it to the controlled diagnostic\.

[Figure4](https://arxiv.org/html/2606.07561#A3.F4)validates the controlled diagnostic against sequential acquisition on three test functions with fitted GP hyperparameters\. The top row shows histograms ofd∂d\_\{\\partial\}for the first point selected after the initial design, aggregated over 1000 independent seeds\. The bottom row shows the controlled diagnostic at the median hyperparameters fitted in the corresponding sequential experiments\. The close agreement between the two rows confirms that the spatial biases visible in the diagnostic are not artifacts of the simplified setup; they persist in full sequential acquisition with MAP\-fitted hyperparameters on real test functions\. Note that the figure labels use “EIG” rather than “VM” because EIG is a monotone transform of the posterior variance; thus, EIG and VM share the same argmax, and their selection profiles are identical\.

![Refer to caption](https://arxiv.org/html/2606.07561v1/x4.png)Figure 4\.Top row:Histograms of boundary distanced∂d\_\{\\partial\}of points selected in sequential acquisition on test functions \(PeaksD=2D\{=\}2, Smooth\-boxD=3D\{=\}3, HartmannD=4D\{=\}4; 1000 seeds each\)\. The dashed line is the geometric density ofd∂d\_\{\\partial\}under uniform sampling from\[0,1\]D\[0,1\]^\{D\}\.Bottom row:Selection profiles from the controlled diagnostic with hyperparameters matched to the sequential experiments \(1000 seeds\)\. Faint dots show raw bin probabilities; solid lines are smoothed; shaded bands are 95% bootstrap confidence intervals\. The dashed baseline is the geometric uniform\-candidate distribution\. VM concentrates at the domain boundary; EPIG concentrates at a characteristic interior shell\.Two\-row figure comparing sequential acquisition \(top\) with the controlled diagnostic \(bottom\) on three test functions across $D=2,3,4$\. Both rows show VM concentrated at the boundary and EPIG concentrated at an interior shell, in close agreement\.
### C\.2\.Lengthscale Robustness

[Figure5](https://arxiv.org/html/2606.07561#A3.F5)extends the controlled diagnostic across a sweep of lengthscalesℓ∈\{0\.20,0\.30,0\.40,0\.55\}\\ell\\in\\\{0\.20,0\.30,0\.40,0\.55\\\}atN=50N=50training points, forD∈\{2,3,4,6\}D\\in\\\{2,3,4,6\\\}\(1000 random training sets per setting\)\. The figure isolates the dependence of each acquisition’s spatial preference on the kernel lengthscale\.

![Refer to caption](https://arxiv.org/html/2606.07561v1/x5.png)Figure 5\.Selection profiles across a lengthscale sweep atN=50N=50training points \(1000 seeds per setting\)\. Rows: dimensionD∈\{2,3,4,6\}D\\in\\\{2,3,4,6\\\}; columns: lengthscaleℓ∈\{0\.20,0\.30,0\.40,0\.55\}\\ell\\in\\\{0\.20,0\.30,0\.40,0\.55\\\}\. Each panel shows the volume\-reweighted acquisition\-argmax distribution for variance maximization \(blue\), negative integrated posterior variance \(red\), and expected predictive information gain \(orange\), against the geometric no\-preference baseline \(dotted gray\)\. VM concentrates at the boundary in every\(D,ℓ\)\(D,\\ell\)setting: the pattern persists across the entire lengthscale range\. NIPV and EPIG sit at the interior in most settings; their boundary fraction grows at largeℓ\\ellin lowDD, where the kernel range becomes a substantial fraction of the domain, and the boundary truncation begins to dominate the global integration\.Four\-by\-four grid of selection\-probability curves across dimensions \(rows: $D=2,3,4,6$\) and lengthscales \(columns: $\\ell=0\.20,0\.30,0\.40,0\.55$\)\. VM concentrates at the boundary in every panel; NIPV and EPIG sit at the interior in most panels, drifting toward the boundary only at the largest $\\ell$ in low $D$\.[Figure6](https://arxiv.org/html/2606.07561#A3.F6)shows normalized mean acquisition\-value profiles across the same sweep: rows areD∈\{2,3,4,6\}D\\in\\\{2,3,4,6\\\}, columns areℓ∈\{0\.20,0\.30,0\.40,0\.55\}\\ell\\in\\\{0\.20,0\.30,0\.40,0\.55\\\}\. In every panel, the normalized mean acquisition value stays close to one acrossd∂d\_\{\\partial\}\(the dotted reference\), even where the corresponding selection profile \([Figure5](https://arxiv.org/html/2606.07561#A3.F5)\) shows a sharp spatial preference, reinforcing the argmax\-versus\-mean point of[SectionB\.1](https://arxiv.org/html/2606.07561#A2.SS1)across the full robustness sweep\.

![Refer to caption](https://arxiv.org/html/2606.07561v1/x6.png)Figure 6\.Normalized mean acquisition\-value profiles across the same sweep as[Figure5](https://arxiv.org/html/2606.07561#A3.F5)\(1000 seeds per setting\)\. Curves are per\-seed normalized so the integral overd∂d\_\{\\partial\}equals0\.50\.5, so a perfectly flat profile sits at1\.01\.0\(the dotted reference\)\. Shaded bands are 95% confidence intervals on the across\-seed mean\. Y\-axes are scaled per panel and centered on1\.01\.0so the within\-panel deviation structure is visible at each panel’s natural scale\. Mean profiles remain narrowly clustered around11in every panel, while the corresponding selection profiles \([Figure5](https://arxiv.org/html/2606.07561#A3.F5)\) show sharp boundary or interior preferences\.Four\-by\-four grid of normalized mean acquisition value versus boundary distance, across dimensions and lengthscales\. Curves remain narrowly clustered around one in every panel, in contrast to the sharp spatial preferences visible in the corresponding selection\-profile figure\.

## Appendix DOne\-Step BO Placement

The selection\-profile diagnostic isolates the observation\-independent family \(VM, NIPV, EPIG\), whose spatial preference, under fixed hyperparameters, is determined by the kernel geometry and the training locations alone \([Remark2\.1](https://arxiv.org/html/2606.07561#S2.Thmremark1)\)\. Full BO acquisitions \(EI, UCB, PI\) add a posterior\-mean component, and the boundary inflation enters only through the variance term\. Whether the bias survives in practice depends on the regime: at smallNNthe mean is uninformative, and exploration dominates; as data accumulates, the mean takes over\.

We extend the diagnostic to a one\-step BO setting on the 2D Gramacy–Lee function\(Surjanovic and Bingham,[2013](https://arxiv.org/html/2606.07561#bib.bib12)\)f​\(x1,x2\)=x1​exp⁡\(−x12−x22\)f\(x\_\{1\},x\_\{2\}\)=x\_\{1\}\\exp\(\-x\_\{1\}^\{2\}\-x\_\{2\}^\{2\}\)\(mapped from\[0,1\]2\[0,1\]^\{2\}to\[−2,6\]2\[\-2,6\]^\{2\}\), which has localized structure in the lower\-left quadrant and is essentially flat elsewhere: a regime where the posterior mean is informative only in a small region\. For eachN0∈\{5,10,25\}N\_\{0\}\\in\\\{5,10,25\\\}and10 00010\\,000Sobol initial designs, we condition a GP with an isotropic Matérn\-5/2\\nicefrac\{\{5\}\}\{\{2\}\}kernel \(fixedℓ=0\.25\\ell=0\.25, fixed noise at5%5\\%ofstd​\(f\)\\mathrm\{std\}\(f\)\) and record the argmax of two acquisitions on a100×100100\\times 100candidate grid: pure varianceVM​\(𝒙\)=s2​\(𝒙\)\\mathrm\{VM\}\(\{\\bm\{x\}\}\)=s^\{2\}\(\{\\bm\{x\}\}\)as the mean\-free reference, andUCB​\(𝒙\)=μ​\(𝒙\)\+κ​s​\(𝒙\)\\mathrm\{UCB\}\(\{\\bm\{x\}\}\)=\\mu\(\{\\bm\{x\}\}\)\+\\kappa\\,s\(\{\\bm\{x\}\}\)withκ=2\\kappa=2\(Srinivaset al\.,[2010](https://arxiv.org/html/2606.07561#bib.bib16)\)\. Sobol initial designs control for the interior\-gap confound: with a well\-spread initial set, there are no large random interior gaps, so a boundary selection reflects the inflation rather than a sparsity accident\.

[Figure7](https://arxiv.org/html/2606.07561#A4.F7)shows the result\. VM places all selections withind∂<0\.05d\_\{\\partial\}<0\.05of the boundary at everyN0N\_\{0\}, the data\-blind corner concentration isolated by the function\-free diagnostic, now confirmed in the presence of training data and a fitted GP\. UCB transitions withN0N\_\{0\}:48%48\\%of selections are near\-boundary atN0=5N\_\{0\}=5,40%40\\%atN0=10N\_\{0\}=10, and15%15\\%atN0=25N\_\{0\}=25, compared with a geometric baseline of∼16%\\sim 16\\%under uniform sampling\. The boundary pull is strongest at smallN0N\_\{0\}, where the posterior mean is weakly informative and the boundary\-inflated variance dominates\. AtN0=25N\_\{0\}=25, UCB is close to the geometric baseline, indicating that the mean component can eventually overcome the variance artifact when enough data are available\.

![Refer to caption](https://arxiv.org/html/2606.07561v1/x7.png)Figure 7\.One\-step BO placement on the 2D Gramacy–Lee function with an isotropic Matérn\-5/2\\nicefrac\{\{5\}\}\{\{2\}\}kernel \(ℓ=0\.25\\ell=0\.25\) and Sobol initial designs of sizeN0N\_\{0\}\.Top row\(VM, posterior variance only\): selection density locked at the corners across allN0N\_\{0\}, independent of the data\.Bottom row\(UCB,μ\+2​σ\\mu\+2\\sigma\): at smallN0N\_\{0\}the mean is uninformative, and the boundary\-inflated variance dominates, with∼48%\\sim 48\\%of selections near the boundary atN0=5N\_\{0\}=5; asN0N\_\{0\}grows, the mean takes over and selection concentrates on the function’s interior structure\. Per\-panel\-normalized densities over10 00010\\,000random initial designs, power\-law color scale\. Function contours shown for reference \(solid = positive levels, dashed = negative\)\.Two\-by\-three grid of 2D selection\-density heatmaps overlaid with function contours\. The top row \(VM\) shows density locked at the four corners for all three initial design sizes\. Bottom row \(UCB\) shows density at the corners at small $N\_\{0\}$, transitioning to the function’s interior structure at larger $N\_\{0\}$\.
## Appendix ENeumann Kernel Correction

The root cause of boundary variance inflation is the truncation of the kernel’s correlation neighborhood at the domain boundary\. A Neumann, or mirror\-image, kernel addresses this by extending the covariance through reflected images, producing an even extension that satisfies Neumann \(zero normal derivative\) boundary conditions on\[0,1\]D\[0,1\]^\{D\}\.

For a stationary one\-dimensional base kernelκ​\(r\)\\kappa\(r\), the Neumann kernel is approximated by the finite image sum

\(7\)kN​\(x,x′\)=∑n=−LL\[κ​\(x−x′\+2​n\)\+κ​\(x\+x′\+2​n\)\],k\_\{N\}\(x,x^\{\\prime\}\)=\\sum\_\{n=\-L\}^\{L\}\\bigl\[\\kappa\(x\-x^\{\\prime\}\+2n\)\+\\kappa\(x\+x^\{\\prime\}\+2n\)\\bigr\],whereLLcontrols the number of image layers\. The first term accounts for periodic translations and the second for reflections about the boundary\. In all experiments,LLis chosen large enough that additional image layers have negligible contribution at the lengthscales considered\. Together, the image terms make the effective correlation neighborhood more symmetric nearx=0x=0andx=1x=1\.

For a product kernel on\[0,1\]D\[0,1\]^\{D\}, we apply this construction dimension\-wise and normalize to unit diagonal:

\(8\)k~N​\(𝒙,𝒙′\)=kN​\(𝒙,𝒙′\)kN​\(𝒙,𝒙\)​kN​\(𝒙′,𝒙′\)\.\\tilde\{k\}\_\{N\}\(\{\\bm\{x\}\},\{\\bm\{x\}\}^\{\\prime\}\)=\\frac\{k\_\{N\}\(\{\\bm\{x\}\},\{\\bm\{x\}\}^\{\\prime\}\)\}\{\\sqrt\{k\_\{N\}\(\{\\bm\{x\}\},\{\\bm\{x\}\}\)\\,k\_\{N\}\(\{\\bm\{x\}\}^\{\\prime\},\{\\bm\{x\}\}^\{\\prime\}\)\}\}\.The normalization is important: without it, the unnormalized kernel has a position\-dependent diagonal that would inflate the prior variance near boundaries, counteracting the correction\. With normalization, the correction arises from off\-diagonal correlations: a boundary point gains correlations with training points on the “reflected” side, increasing the effective number of contributing neighbors\.

This is not a claim that the unknown function satisfies Neumann boundary conditions\. Rather, the mirror\-image construction corrects the geometric artifact in the posterior variance that causes boundary oversampling\.

[Figure8](https://arxiv.org/html/2606.07561#A5.F8)shows the diagnostic’s response to swapping the standard product Matérn kernel for the Neumann kernel, at the same controlled\-diagnostic settings as[Figure2](https://arxiv.org/html/2606.07561#S4.F2)\(N=50N=50,ℓ\\ellperDDas listed\)\. AcrossD∈\{2,3,4,6\}D\\in\\\{2,3,4,6\\\}, the Neumann kernel reduces the fraction of VM selections withind∂<0\.05d\_\{\\partial\}<0\.05of the boundary by 6 to 15 percentage points, with the reduction growing withDD\. Expressed as a fraction of the excess over the geometric no\-preference baseline, the correction removes7%7\\%atD=2D=2,12%12\\%atD=3D=3,15%15\\%atD=4D=4, and28%28\\%atD=6D=6, consistent with the2D2^\{D\}corner amplification of the product construction\. The baseline shown is the geometricP​\(d∂​<0\.05∣​d∂∈\[0\.01,0\.5\]\)P\(d\_\{\\partial\}<0\.05\\mid d\_\{\\partial\}\\in\[0\.01,0\.5\]\)to match the candidate grid’s range; selections use the volume\-reweighted argmax of[Figure2](https://arxiv.org/html/2606.07561#S4.F2)\.

![Refer to caption](https://arxiv.org/html/2606.07561v1/x8.png)Figure 8\.Selection\-profile response to the Neumann correction at the matched controlled\-diagnostic settings of[Figure2](https://arxiv.org/html/2606.07561#S4.F2)\(N=50N=50\)\. Fraction of acquisition argmax withind∂<0\.05d\_\{\\partial\}<0\.05under: the geometric no\-preference baseline \(gray\); VM with the standard product Matérn kernel \(blue\); VM with the Neumann kernel \(green\); and NIPV \(red\) and EPIG \(orange\) with the standard kernel for context\. VM has a strong boundary bias, and Neumann shifts it toward the baseline, with the reduction increasing asDDincreases\. NIPV is intermediate atD=2D=2and self\-corrects at higherDD; EPIG sits below the baseline throughout\. Both penalize boundary candidates via their respective integration structures, so the kernel modification is most relevant to VM\.Bar chart with four groups of five bars showing percent near\-boundary selections at dimensions 2, 3, 4, and 6 for: a no\-preference geometric baseline, VM with the standard kernel, VM with the Neumann kernel, NIPV, and EPIG\. VM stays at 94\-99% across all dimensions; the Neumann variant shifts it downward by 5\-15 percentage points\.For context,[Figure8](https://arxiv.org/html/2606.07561#A5.F8)also shows NIPV and EPIG at the same settings, with the standard kernel\. EPIG sits at∼0%\\sim 0\\%in every dimension, well below the baseline: its averaging over targets already penalizes boundary candidates\. NIPV is intermediate atD=2D=2\(15%15\\%, just below the16%16\\%baseline\) and drops to∼0%\\sim 0\\%forD∈\{3,4,6\}D\\in\\\{3,4,6\\\}, matching the “partial correction through global integration” picture from[SectionA\.3](https://arxiv.org/html/2606.07561#A1.SS3)\. The kernel\-level Neumann modification is therefore most consequential for VM, the acquisition whose argmax most directly inherits the boundary inflation\.

The correction is partial in every dimension\. Reflection modifies the covariance geometry but does not introduce independent observations outside the domain, so residual boundary preference persists\. Fully equalizing the variance field may require nonstationary kernels that explicitly encode the boundary, or input warping\(Snoeket al\.,[2014](https://arxiv.org/html/2606.07561#bib.bib15)\)that maps boundary regions inward; both directions are left for follow\-up work\. Here, the Neumann result serves primarily as a causal test of the diagnosis: modifying the truncated correlation neighborhood shifts the selection profile in the predicted direction, confirming that the bias is geometric in origin\.