#  Artificial Intelligence and Human Geography

Chapter in the Encyclopedia of Human Geography

Song Gao<EMAIL_ADDRESS>Geospatial Data Science Lab, Department of
Geography, University of Wisconsin-Madison, 550 N Park Stret, Madison, 53706,
WI, USA

###### Abstract

This paper examines the recent advances and applications of AI in human
geography especially the use of machine (deep) learning, including place
representation and modeling, spatial analysis and predictive mapping, and
urban planning and design. AI technologies have enabled deeper insights into
complex human-environment interactions, contributing to more effective
scientific exploration, understanding of social dynamics, and spatial
decision-making. Furthermore, human geography offers crucial contributions to
AI, particularly in context-aware model development, human-centered design,
biases and ethical considerations, and data privacy. The synergy beween AI and
human geography is essential for addressing global challenges like disaster
resilience, poverty, and equitable resource access. This interdisciplinary
collaboration between AI and geography will help advance the development of
GeoAI and promise a better and sustainable world for all.

###### keywords:

AI, GeoAI, human-in-the-loop, ethics, trustworthiness

##  1 Introduction

Artificial intelligence (AI) is a field in computer science and engineering
that focuses on developing intelligent machines capable of performing problem-
solving tasks and achieving goals that typically require human intelligence
[64, 45]. There are generally two types of AI: Weak AI and Strong AI. Weak AI
sits at the foundation of AI development, with its focused human-like specific
abilities (e.g., conversation in Siri) and practical applications (e.g.,
facial recognition from images and videos) in specific contexts, while Strong
AI (known as Artificial General Intelligence) represents the ultimate goal of
AI research, aiming to replicate or even surpass human intelligence to solve
problems, learn, and plan for the future with self-aware consciousness [18].

AI is revolutionizing many domains including geographical sciences, presenting
both immense opportunities and formidable challenges. The rapid advancement of
AI is fueled by theoretical breakthroughs in neural science and computer
science, the ubiquity of big data, advanced computer hardware such as graphics
processing units (GPUs) and powerful high-performance computing platforms that
enable the efficient development training, and deployment of AI models, and
the flourishing AI products, services and applications that are changing human
behaviors [9, 39, 62].

As [5] expressed her excitement about the AI-driven evolution—“It is not often
that geography is touched by a development having the potential to affect
substantially all of the practical, technical, methodological, theoretical and
philosophical aspects of our work.” The intersection of AI and geography has a
rich and evolving history, with its early roots documented in works by [5, 60,
51]. Prior to the surge of deep learning research in the 2010s [38] and its
applications in geography and earth sciences [24, 19, 40, 55], significant AI
advancements included theoretical explorations in the 1950s and 1960s [3]; the
emergence of artificial neural networks, heuristic search algorithms,
knowledge-based symbolic expert systems, neurocomputing, and artificial life
concepts (e.g., cellular automata) in the 1980s; the development of genetic
programming, fuzzy logics, and hybrid intelligent systems in the 1990s [51];
the integration of ontology and semantics for information retrieval and
knowledge graphs in the 2000s; and the recent development of foundation model
[1, 44]. All of these advances have laid the foundation for the burgeoning
field of Geospatial artificial intelligence (GeoAI) [15].

Rooted in geography and geographic information science (GIScience), GeoAI has
now emerged as a transformative interdisciplinary field, integrating
geographical studies with AI techniques, particularly spatially-explicit
machine learning and deep learning methods as well as geo-knowledge graphs [9,
28, 43]. This convergence has spurred groundbreaking advancements in both
academia and industry. GeoAI encompasses the development of intelligent
computer programs that emulate human perception, spatial reasoning, and the
ability to uncover insights from geographical phenomena and understanding
their complex dynamics through coupling spatial knowledge and earth science
process-based models with deep neural networks [10, 55]. The researches in
GeoAI can deepen our understanding of complex human-environmental systems and
their interactions, particularly within a spatial context.

Regarding the topics of AI in human geography, a recent conversation about
GeoAI, counter-AI, and human geography by [29] discussed opportunities and
challenges in defining intelligence for machines, and the role of humans in AI
development and decision-making processes, emphasizing the need for ethical
considerations and accountability in the development of GeoAI. In the
following, I will first introduce some of the recent advances of AI in human
geography and then discuss the role of human geography in AI development.

##  2 AI in Human Geography

Human geography has a long tradition of studying the spatially and temporally
varying interrelationships between people, place, and environment, and the
ways in which people and their activities shape, and are shaped by, the
physical, economic, political, and cultural environment [30]. The integration
of AI technologies has significantly expanded the capacity of researchers and
practitioners in understanding the complex human-environment interactions.
Geographers have explored a wide range of topics and made advances in the
following areas (but not limited to).

Place Representation and Modeling:  Place can serve as a function between
location and people [47], a function of location, activity and time [46], and
a function of social relations [17]. Traditionally, data about places were
collected through mapping agencies (e.g., digital gazetteers) and survey-based
narratives. The emergence of geospatial big data brings new opportunities to
extract fine spatiotemporal resolution of human-place interaction data and
understand the rich place semantics from large-scale volunteered geographic
information and crowdsourced data streams, such as social media posts
(including texts, photos, and videos), GPS tracks, text reviews on points of
interest (POIs) and neighborhoods, and other Web documents [12, 13, 41]. A
wide-range of AI methods provide new opportunities to understand and extract
the characteristics of places as well as associated human activities,
experiences, emotions, and movements in different contexts. For example, urban
areas of interest that attract people’s attention but with different spatial
patterns and region-specific semantics were automatically extracted from
georeferenced photos posted on social media by employing spatial and spectral
clustering algorithms, and natural language processing techniques [22]. In
addition, the spatial and hierarchical semantics between places that support
human cognition of places, identification of urban functional regions, and
information retrieval in digital maps, were modelled using POI embeddings
(i.e., multidimensional vectors) via topic modeling [11] and deep learning
techniques [26]. Patterns and relations between places can be computed and
extracted from collective human descriptions using place graphs [4]. Empowered
by the state-of-the-art human face and emotion recognition AI techniques,
human emotions at different places were extracted from millions of
georeferenced photos [31], which enabled the study of spatially embodied
emotions in human geography [59]. Recently, researchers found that geo-
knowledge-guided large language models (e.g., ChatGPT) improved the extraction
of location textual descriptions from disaster-related social media messages,
which can facilitate collaborations between disaster response experts and AI
developers and ultimately help save human lives [25].

Spatial Analysis and Predictive Mapping:  AI has revolutionized geospatial
data analysis, interpretation, and modeling of human-environment relations by
providing innovative spatial analytical tools and methodologies. New or
modification of spatial analysis methods and spatial statistical models have
been proposed by incorporating deep neural networks, which we termed as
Intelligent Spatial Analytics [73]. Geographically weighted artificial neural
network (GWANN) [20] and geographically neural network weighted regression
(GWNNR) [7] models were developed by constructing the spatial nonstationary
weights with neural networks and feeding into the geographically weighted
regression for more accurate predictive modeling. Similarly, new spatial
regression graph convolutional neural networks (SRGCNNs) were developed for
the spatial regression analysis of geographical multivariate distributions and
outperformed traditional spatial autoregressive models [74]. Machine learning
and deep learning algorithms have also been widely used for data-driven
modeling and predictive mapping in geography [48] and identifying patterns and
relationships that would be nearly impossible for human researchers to discern
manually [37]. For example, [72] applied the graph convolutional neural
networks (GCNNs) for modeling place graphs and inferring the unknown
properties of a place with a high prediction accuracy by utilizing both the
observed attributes (e.g., visual and functional features) and the relational
characteristics (e.g., distance, adjacency, spatial interaction intensity) of
the places. When applying state-of-the-art deep learning models for population
mapping in geography, the issues on model selection, neighboring effects, and
systematic biases play an important role on the predictive mapping accuracy
[27]. Researchers found that population estimation performance reduced with
increasing neighborhood sizes and a pervasive bias existed in which all tested
deep learning models (i.e., VGG, ResNet, Xception, and DenseNet) overestimated
sparsely populated patches and underestimated densely populated ones [27].
This work also demonstrated the importance of spatial concepts and some
limitations when applying AI methods in human geography due to spatial
heterogeneity. Future efforts need to be made towards the (weak) replication
of GeoAI models across space and time in the social and environmental sciences
[19].

Urban Planning and Design:  Human geographers have long been interested in
understanding the intricate relationships between people and the environment
and collaborating with urban planners for planning our cities. By applying AI
algorithms, they can effectively interpret large geospatial datasets,
including census surveys, satellite imagery, street-view images, social media
posts, and other sensor data to support urban analytics [6]. By utilizing the
state-of-the-art deep learning models, human perception of neighborhood
playability for childhood development in cities [36] and urban visual
intelligence about the hidden neighborhood socioeconomic status, such as
poverty status and health outcomes, can be extracted from street-view images
[8, 32]; near real-time global land cover and land use patterns and temporal
changes can be detected from 10m spatial resolution of Sentinel-2 remote
sensing imagery [2]. Urban geographers and planners can also extract valuable
insights into population dynamics, traffic patterns, land use patterns,
environmental changes, and the interactions between social and ecological
systems, so as to optimize urban spatial structure, improving transportation
efficiency and quality of living, and environmental sustainability [50, 40].
AI-driven predictive modeling can assist in rooftop solar potential estimation
for sustainable city design [66], and anticipating changes in population
density and mobility, facilitating better resource allocation and disaster
preparedness in [75, 65]. Additionally, AI-powered semantic and sentiment
analyses of online neighborhood reviews helped better understand the
perceptions and emotions of people toward their living neighborhoods and
environments, which can be applied for supporting urban planning and improving
quality of life studies [23]. Furthermore, the advances in generative AI
models brings exiting opportunities for urban planners to facilitate the
automatic rendering of urban master plans [68] and building floorplan layouts
[67] via generative adversarial networks (GANs), otherwise they need to rely
on subjective design and labor-intensive production process. However, the
generative AI also introduces the deep fakes into geography [71] and we need
to make more efforts on the investigation of trustworthiness of such AI
generated geospatial data and city plans.

##  3 Human Geography for AI

Human geography also contributes to the development and applications of AI by
providing valuable insights into the social, cultural, and economic factors
that shape human behavior and interactions with the technology.

Context-Aware and Debiasing AI:  Human geography contributes to AI through the
development of context-aware models and the analysis of biases of AI models
when applied in heterogeneous contexts and places. Human perception of the
physical and socioeconomic environments through deep learning has attracted
large interest of scholars from geography, computer sciences, urban planning,
and environmental psychology [16, 8, 32, 21]. Human perception of places
reflects people’s subjective feelings towards their surrounding environment,
e.g., whether a place is safe, lively, and wealthy [70]. However, there exist
multiple types of biases in existing AI-driven computational frameworks when
quantifying human perceptions of places: data bias, model bias, and perception
bias. [33] compared the safety perception measures from the survey based on
neighborhood residents’ responses with those from the GeoAI approach with
street-view imagery to better understand the relationship between the two
types of measures. This research found that citywide residents, but not
neighborhood residents, may feel economically vibrant places look safe; elder
people may underestimate the safety of their living places which may enlarge
perception bias. Therefore, integrating cultural, historical, and
socioeconomic insights into pre-trained AI models and fine-tuning them with
local data and participant inputs may enhance AI’s ability to comprehend and
respond to spatial heterogeneity and diverse human experiences. This approach
not only improves the accuracy of AI applications but also fosters inclusivity
and respect for social diversity across geographic regions.

Data Privacy and Ethical Framework:  Human geography addresses ethical
concerns related to data collection and usage, as well as community trust
[57]. This knowledge is invaluable for AI developers who must navigate the
ethical implications of collecting and analyzing vast amounts of data
regarding privacy and security. In GeoAI development, such concerns often
refer to the use or exposure of sensitive geospatial information such as a
user’s home location, workspace, POI preferences, daily trajectories, and
associated personal inferences based on such information [35, 54]. There is a
trade-off between the data utility in downstream tasks (e.g., user profiling)
and user privacy [14]. In human mobility studies, common user privacy
protection methods include de-identification, geomasking, trajectory
k-anonymization, and differential privacy. Recently, deep neural networks such
as LSTM and GANs have been used in human movement trajectory generation tasks
with the aim to balance the data utility-privacy trade-off [52, 54]. With the
substantial advances in AI foundation models, recent studies, however,
revealed that the development and use of foundation models could potentially
unveil substantial privacy and security risks, including the disclosure of
sensitive information, representational bias, hallucinations, and misuse [1].
In the lifecycle of building and using GeoAI foundation models [44], a series
of potential privacy and security risks that existed around the pre-training
and fine-tuning stages with geospatial data, centralized serving and tooling,
prompting-based interaction, and feedback mechanisms [53]. Looking ahead,
human geographers and GIScience scholars can further contribute to
establishing guidelines that ensure responsible and privacy-conscious GeoAI
practices.

Human-Centered Design:  Understanding the perspectives, experiences, emotions,
and lived realities of people in different places and regions is essential for
human geography research [63]. Therefore, human geography emphasizes the
importance of human-centered design in place-based innovation and policy
making, aligning seamlessly with the principles of ethical (Geo-)AI
development [58, 34] and user-centered design in cartography [56]. By
incorporating human geographers in the (Geo-)AI design processes, developers
can create AI systems that better serve diverse communities and mitigate
potential negative impacts on our society, which is well aligned with the core
idea of keeping “human-in-the-loop” for the development of AI systems and
algorithms [29, 69]. Several principles and practices like collaborative
engagement with data [42], accountability mechanisms in AI systems design [49]
and empathy [61] can be considered for addressing the social-technical
challenges in the AI paradigm.

##  4 Conclusion

The relationship between AI and human geography is not a one-way street. As AI
becomes increasingly sophisticated, its applications in human geography will
become even more diverse and impactful. At the same time, the abovementioned
critical insights provided by human geographers will be crucial for ensuring
that less-biased (Geo-)AI is developed and deployed in a responsible and
ethical manner. Mitigating multiple types of biases throughout the “data-
model-action” loop in human-environment systems may ensure equitable outcomes
and AI-driven spatial decision making. By working together, AI and human
geography can address some of the most pressing challenges facing human-
environment nexus today, including climate change, poverty, food security,
health, disaster resilience, and equitable access to resources and services.
This interdisciplinary collaboration has the potential to create a better and
sustainable world for all.

## References

  *   * Bommasani et al [2021] Bommasani R, Hudson DA, Adeli E, et al (2021) On the opportunities and risks of foundation models. arXiv preprint arXiv:210807258
  * Brown et al [2022] Brown CF, Brumby SP, Guzder-Williams B, et al (2022) Dynamic world, near real-time global 10 m land use land cover mapping. Scientific Data 9(1):251
  * Buchanan [2005] Buchanan BG (2005) A (very) brief history of artificial intelligence. AI Magazine 26(4):53–53
  * Chen et al [2018] Chen H, Winter S, Vasardani M (2018) Georeferencing places from collective human descriptions using place graphs. Journal of Spatial Information Science (17):31–62
  * Couclelis [1986] Couclelis H (1986) Artificial intelligence in geography: Conjectures on the shape of things to come. The Professional Geographer 38(1):1–11
  * De Sabbata et al [2023] De Sabbata S, Ballatore A, Miller HJ, et al (2023) Geoai in urban analytics
  * Du et al [2020] Du Z, Wang Z, Wu S, et al (2020) Geographically neural network weighted regression for the accurate estimation of spatial non-stationarity. International Journal of Geographical Information Science 34(7):1353–1377
  * Fan et al [2023] Fan Z, Zhang F, Loo BP, et al (2023) Urban visual intelligence: Uncovering hidden city profiles with street view images. Proceedings of the National Academy of Sciences 120(27):e2220417120
  * Gao [2020] Gao S (2020) A review of recent researches and reflections on geospatial artificial intelligence. Geomatics and Information Science of Wuhan University 45(12):1865–1874
  * Gao [2021] Gao S (2021) Geospatial artificial intelligence (GeoAI). Oxford University Press
  * Gao et al [2017a] Gao S, Janowicz K, Couclelis H (2017a) Extracting urban functional regions from points of interest and human activities on location-based social networks. Transactions in GIS 21(3):446–467
  * Gao et al [2017b] Gao S, Janowicz K, Montello DR, et al (2017b) A data-synthesis-driven method for detecting and extracting vague cognitive regions. International Journal of Geographical Information Science 31(6):1245–1271
  * Gao et al [2017c] Gao S, Li L, Li W, et al (2017c) Constructing gazetteers from volunteered big geo-data based on hadoop. Computers, Environment and Urban Systems 61:172–186
  * Gao et al [2019] Gao S, Rao J, Liu X, et al (2019) Exploring the effectiveness of geomasking techniques for protecting the geoprivacy of twitter users. Journal of Spatial Information Science (19):105–129
  * Gao et al [2023] Gao S, Hu Y, Li W (2023) Handbook of Geospatial Artificial Intelligence. CRC Press, URL https://doi.org/10.1201/9781003308423
  * Gebru et al [2017] Gebru T, Krause J, Wang Y, et al (2017) Using deep learning and google street view to estimate the demographic makeup of neighborhoods across the united states. Proceedings of the National Academy of Sciences 114(50):13108–13113
  * Giordano and Cole [2018] Giordano A, Cole T (2018) The limits of gis: Towards a gis of place. Transactions in GIS 22(3):664–676
  * Goertzel [2014] Goertzel B (2014) Artificial general intelligence: concept, state of the art, and future prospects. Journal of Artificial Intelligence 5(1):1
  * Goodchild and Li [2021] Goodchild MF, Li W (2021) Replication across space and time must be weak in the social and environmental sciences. Proceedings of the National Academy of Sciences 118(35):e2015759118
  * Hagenauer and Helbich [2022] Hagenauer J, Helbich M (2022) A geographically weighted artificial neural network. International Journal of Geographical Information Science 36(2):215–235
  * Han et al [2022] Han X, Wang L, Seo SH, et al (2022) Measuring perceived psychological stress in urban built environments using google street view and deep learning. Frontiers in Public Health 10:891736
  * Hu et al [2015] Hu Y, Gao S, Janowicz K, et al (2015) Extracting and understanding urban areas of interest using geotagged photos. Computers, Environment and Urban Systems 54:240–254
  * Hu et al [2019a] Hu Y, Deng C, Zhou Z (2019a) A semantic and sentiment analysis on online neighborhood reviews for understanding the perceptions of people toward their living environments. Annals of the American Association of Geographers 109(4):1052–1073
  * Hu et al [2019b] Hu Y, Li W, Wright D, et al (2019b) Artificial intelligence approaches. The Geographic Information Science & Technology Body of Knowledge
  * Hu et al [2023] Hu Y, Mai G, Cundy C, et al (2023) Geo-knowledge-guided gpt models improve the extraction of location descriptions from disaster-related social media messages. International Journal of Geographical Information Science 37(11):2289–2318
  * Huang et al [2022] Huang W, Cui L, Chen M, et al (2022) Estimating urban functional distributions with semantics preserved poi embedding. International Journal of Geographical Information Science 36(10):1905–1930
  * Huang et al [2021] Huang X, Zhu D, Zhang F, et al (2021) Sensing population distribution from satellite imagery via deep learning: Model selection, neighboring effects, and systematic biases. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 14:5137–5151
  * Janowicz et al [2020] Janowicz K, Gao S, McKenzie G, et al (2020) GeoAI: spatially explicit artificial intelligence techniques for geographic knowledge discovery and beyond. International Journal of Geographical Information Science 34(4):625–636
  * Janowicz et al [2022] Janowicz K, Sieber R, Crampton J (2022) Geoai, counter-ai, and human geography: A conversation. Dialogues in Human Geography 12(3):446–458
  * Jones [2012] Jones A (2012) Human geography: The basics. Routledge
  * Kang et al [2019] Kang Y, Jia Q, Gao S, et al (2019) Extracting human emotions at different places based on facial expressions and spatial clustering analysis. Transactions in GIS 23(3):450–480
  * Kang et al [2020] Kang Y, Zhang F, Gao S, et al (2020) A review of urban physical environment sensing using street view imagery in public health studies. Annals of GIS 26(3):263–275
  * Kang et al [2023] Kang Y, Abraham J, Ceccato V, et al (2023) Assessing differences in safety perceptions using geoai and survey across neighbourhoods in stockholm, sweden. Landscape and Urban Planning 236:104768
  * Kang et al [2024] Kang Y, Gao S, Roth R (2024) Artificial intelligence studies in cartography: A review and synthesis of methods, applications, and ethics. pp 1–20
  * Keßler and McKenzie [2018] Keßler C, McKenzie G (2018) A geoprivacy manifesto. Transactions in GIS 22(1):3–19
  * Kruse et al [2021] Kruse J, Kang Y, Liu YN, et al (2021) Places for play: Understanding human perception of playability in cities using street view images and deep learning. Computers, Environment and Urban Systems 90:101693
  * Lavallin and Downs [2021] Lavallin A, Downs JA (2021) Machine learning in geography–past, present, and future. Geography Compass 15(5):e12563
  * LeCun et al [2015] LeCun Y, Bengio Y, Hinton G (2015) Deep learning. nature 521(7553):436–444
  * Li [2020] Li W (2020) GeoAI: Where machine learning and big data converge in GIScience. Journal of Spatial Information Science (20):71–77
  * Liu and Biljecki [2022] Liu P, Biljecki F (2022) A review of spatially-explicit geoai applications in urban geography. International Journal of Applied Earth Observation and Geoinformation 112:102936
  * Liu et al [2015] Liu Y, Liu X, Gao S, et al (2015) Social sensing: A new approach to understanding our socioeconomic environments. Annals of the Association of American Geographers 105(3):512–530
  * Lloyd and Dykes [2011] Lloyd D, Dykes J (2011) Human-centered approaches in geovisualization design: Investigating multiple methods through a long-term case study. IEEE Transactions on Visualization and Computer Graphics 17(12):2498–2507
  * Mai et al [2022] Mai G, Hu Y, Gao S, et al (2022) Symbolic and subsymbolic geoai: Geospatial knowledge graphs and spatially explicit machine learning
  * Mai et al [2023] Mai G, Huang W, Sun J, et al (2023) On the opportunities and challenges of foundation models for geospatial artificial intelligence. arXiv preprint arXiv:230406798 pp 1–32
  * McCarthy [2004] McCarthy J (2004) What is artificial intelligence? Stanford University pp 1–15
  * McKenzie and Adams [2017] McKenzie G, Adams B (2017) Juxtaposing thematic regions derived from spatial and platial user-generated content. In: 13th international conference on spatial information theory (COSIT 2017), Schloss Dagstuhl-Leibniz-Zentrum fuer Informatik
  * Mennis and Mason [2016] Mennis J, Mason MJ (2016) Modeling place as a relationship between a person and a location. In: International Conference on GIScience Short Paper Proceedings
  * Miller and Goodchild [2015] Miller HJ, Goodchild MF (2015) Data-driven geography. GeoJournal 80:449–461
  * Mittelstadt [2019] Mittelstadt B (2019) Principles alone cannot guarantee ethical ai. Nature Machine Intelligence 1(11):501–507
  * Mortaheb and Jankowski [2023] Mortaheb R, Jankowski P (2023) Smart city re-imagined: City planning and geoai in the age of big data. Journal of Urban Management 12(1):4–15
  * Openshaw and Openshaw [1997] Openshaw S, Openshaw C (1997) Artificial intelligence in geography. John Wiley & Sons, Inc.
  * Rao et al [2020] Rao J, Gao S, Kang Y, et al (2020) Lstm-trajgan: A deep learning approach to trajectory privacy protection. In: 11th International Conference on Geographic Information Science (GIScience 2021), Schloss Dagstuhl-Leibniz-Zentrum für Informatik, p 12
  * Rao et al [2023a] Rao J, Gao S, Mai G, et al (2023a) Building privacy-preserving and secure geospatial artificial intelligence foundation models. In: The 31st ACM International Conference on Advances in Geographic Information Systems (SIGSPATIAL’23). ACM, pp 1–4
  * Rao et al [2023b] Rao J, Gao S, Zhu S (2023b) Cats: Conditional adversarial trajectory synthesis for privacy-preserving trajectory data publication using deep learning approaches. International Journal of Geographical Information Science 37(12):2538–2574
  * Reichstein et al [2019] Reichstein M, Camps-Valls G, Stevens B, et al (2019) Deep learning and process understanding for data-driven earth system science. Nature 566(7743):195–204
  * Roth [2017] Roth R (2017) User interface and user experience (ui/ux) design. Geographic Information Science & Technology Body of Knowledge 2017(Q2):1–10
  * Rowe [2021] Rowe F (2021) Big data and human geography. SocArXiv Papers 2021(0):1–10
  * Siau and Wang [2020] Siau K, Wang W (2020) Artificial intelligence (ai) ethics: ethics of ai and ethical ai. Journal of Database Management (JDM) 31(2):74–87
  * Simonsen et al [2007] Simonsen K, et al (2007) Practice, spatiality and embodied emotions: An outline of a geography of practice. Human Affairs (2):168–181
  * Smith [1984] Smith TR (1984) Artificial intelligence and its applicability to geographical problem solving. The Professional Geographer 36(2):147–158
  * Srinivasan and González [2022] Srinivasan R, González BSM (2022) The role of empathy for artificial intelligence accountability. Journal of Responsible Technology 9:100021
  * Torrens [2018] Torrens PM (2018) Artificial intelligence and behavioral geography. Handbook of Behavioral and Cognitive Geography 20(0):357–371
  * Tuan [1979] Tuan YF (1979) Space and place: humanistic perspective. In: Philosophy in geography. Springer, p 387–427
  * Turing [1950] Turing AM (1950) Computing machinery and intelligence. Mind 49(0):433–460
  * Vongkusolkit and Huang [2021] Vongkusolkit J, Huang Q (2021) Situational awareness extraction: a comprehensive review of social media data classification during natural hazards. Annals of GIS 27(1):5–28
  * Wu and Biljecki [2021] Wu AN, Biljecki F (2021) Roofpedia: Automatic mapping of green and solar roofs for an open roofscape registry and evaluation of urban sustainability. Landscape and Urban Planning 214:104167
  * Wu et al [2022] Wu AN, Stouffs R, Biljecki F (2022) Generative adversarial networks in the built environment: A comprehensive review of the application of gans across data types and scales. Building and Environment p 109477
  * Ye et al [2022] Ye X, Du J, Ye Y (2022) Masterplangan: Facilitating the smart rendering of urban master plans via generative adversarial networks. Environment and Planning B: Urban Analytics and City Science 49(3):794–814
  * Zanzotto [2019] Zanzotto FM (2019) Human-in-the-loop artificial intelligence. Journal of Artificial Intelligence Research 64:243–252
  * Zhang et al [2018] Zhang F, Zhou B, Liu L, et al (2018) Measuring human perceptions of a large-scale urban region using machine learning. Landscape and Urban Planning 180:148–160
  * Zhao et al [2021] Zhao B, Zhang S, Xu C, et al (2021) Deep fake geography? when geospatial data encounter artificial intelligence. Cartography and Geographic Information Science 48(4):338–352
  * Zhu et al [2020] Zhu D, Zhang F, Wang S, et al (2020) Understanding place characteristics in geographic contexts through graph convolutional neural networks. Annals of the American Association of Geographers 110(2):408–420
  * Zhu et al [2022a] Zhu D, Gao S, Cao G (2022a) Towards the intelligent era of spatial analysis and modeling. In: Proceedings of the 5th ACM SIGSPATIAL International Workshop on AI for Geographic Knowledge Discovery. ACM, pp 10–13
  * Zhu et al [2022b] Zhu D, Liu Y, Yao X, et al (2022b) Spatial regression graph convolutional neural networks: A deep learning paradigm for spatial multivariate distributions. GeoInformatica 26(0):645––676
  * Zou et al [2022] Zou L, Cai H, Mandal D, et al (2022) Empowering smart disaster response with social media and geoai. In: AGU Fall Meeting, vol 2022. American Geophysical Union, pp NH15C–0329


# Scalable Bayesian Optimization Using Vecchia Approximations of Gaussian
Processes

Firstname1 Lastname1  Firstname2 Lastname2  Firstname3 Lastname3  Firstname4
Lastname4  Firstname5 Lastname5  Firstname6 Lastname6  Firstname7 Lastname7
Firstname8 Lastname8  Firstname8 Lastname8

###### Abstract

Bayesian optimization is a technique for optimizing black-box target
functions. At the core of Bayesian optimization is a surrogate model that
predicts the output of the target function at previously unseen inputs to
facilitate the selection of promising input values. Gaussian processes (GPs)
are commonly used as surrogate models but are known to scale poorly with the
number of observations. We adapt the Vecchia approximation, a popular GP
approximation from spatial statistics, to enable scalable high-dimensional
Bayesian optimization. We develop several improvements and extensions,
including training warped GPs using mini-batch gradient descent, approximate
neighbor search, and selecting multiple input values in parallel. We focus on
the use of our warped Vecchia GP in trust-region Bayesian optimization via
Thompson sampling. On several test functions and on two reinforcement-learning
problems, our methods compared favorably to the state of the art.

Bayesian optimization, Gaussian proecess, Vecchia Approximation, Sparse
Cholesky, Machine Learning

##  1 Introduction

Bayesian optimization (BO) is a class of techniques for black-box function
optimization (Mockus1989). BO has found success in a variety of areas,
including reinforcement learning (e.g., Calandra2016), tuning of machine-
learning algorithms (e.g., Swersky2013), A/B testing (e.g., Letham2019) and
material discovery (e.g., gomez2018automatic).

The basic idea of BO is to approximate the black-box function with a cheap-to-
evaluate surrogate, typically a GP, and perform optimization with the
surrogate. Starting with an initial collection of function evaluations, the
trained GP is combined with a heuristic to generate an acquisition function,
which assigns values to unseen inputs. The input with the highest acquisition
value is evaluated and the result is added to the dataset. This process
continues iteratively until a predefined stopping criterion is met or the
computational budget is depleted. For a more complete review of Bayesian
optimization see Frazier2018AOptimization.

As most BO techniques use a GP surrogate, they inherit the scalability issues
of GPs. It becomes expensive to search a high-dimensional space when
optimizing the acquisition function.

This work addresses this scalability issue by extending the Vecchia GP
approximation from spatial statistics (e.g., Vecchia1988; Katzfuss2017a) to
work in high dimensions with many observations. Through these extensions, we
demonstrate how Vecchia GPs can be used in concert with existing BO
procedures. To be precise, our contributions are:

  * •

Introduce Vecchia GP approximations to the machine-learning community.

  * •

Adapt Vecchia GPs to BO and improve their flexibility and scaling with warped
kernels, vector quantization, approximate nearest neighbors, minibatch
training and variance correction.

  * •

Empirically demonstrate the utility of Vecchia GPs within existing BO
frameworks.

##  2 Existing Solutions

There has been interest in scaling BO in two ways: higher input dimensions and
more function evaluations. While we mostly focus on the latter here, these two
aspects are related, in that successful high-dimensional optimization
typically requires many function evaluations.

High-dimensional BO is thoroughly studied with many proposed solutions,
including transformations to a lower-dimensional space (gomez2018automatic),
assuming an additive form of the function (kandasamy2015high), and using
sparsity-inducing priors (eriksson2021high).

Scaling BO to many function evaluations is less explored. Existing approaches
rely on additivity assumptions (e.g., wang2018batched) or GP approximations
(mcintire2016sparse; maddox2021conditioning). The Trust Region BO method
(TuRBO; eriksson2019scalable) has also seen great success by relying on
Thompson sampling and limiting the input search to trust regions, which are
determined by the lengthscale parameters of the GP.

Here we argue and show that GP-based BO methods can be made scalable via the
Vecchia approximation, and we establish that existing scalable BO methods can
be improved by using Vecchia GPs instead of inducing point methods (e.g.,
Hensman2013GaussianData), which can suffer from important limitations
(Stein2013a). We hope to provide readers with a set of techniques for
incorporating Vecchia GPs into existing BayesOpt procedures to improve their
scalability.

##  3 Background

###  3.1 Review of Bayesian Optimization

Bayesian optimization aims to solve the problem

$\displaystyle\min_{\bm{x}\in\mathcal{X}}f(x),$ (1)

for a “black-box” objective function $f$ (called the target function) over
some input space $\mathcal{X}\subset\mathbb{R}^{d}$. In BayesOpt, $f$ is often
modeled as a Gaussian process (GP), $f(\cdot)\sim
GP(\mu_{\bm{\theta}},K_{\bm{\theta}})$, whose mean $\mu_{\bm{\theta}}$ and
kernel function $K_{\bm{\theta}}$ depend on hyperparameters $\bm{\theta}$.

Function evaluations $y_{i}=f(\bm{x}_{i})$ are obtained sequentially. Given
the first $n$ observations $\bm{y}_{1:n}$ at inputs $\bm{x}_{1:n}$, the goal
is to find a new input $\bm{x}_{n+1}$, such that the global optimum of $f$ is
likely at $\bm{x}_{n+1}$.

This is usually achieved by maximizing an acquisition function,
$\bm{x}_{n+1}=\arg\max_{\bm{x}\in\mathcal{X}}a_{\bm{\theta}^{*}}(\bm{x}|\bm{y}_{1:n})$,
which is a function that assigns a measure of value to unseen locations based
on a utility function, and where $\bm{\theta}^{*}=\arg\max_{\bm{\theta}}\log
p_{\theta}(\bm{y}_{1:n})$ is the maximum likelihood estimate.

Alternatively, new inputs can be chosen using Thompson sampling
(Thompson1933), which has been shown to work well empirically (Chapelle2011).
In Thompson sampling, we consider a collection of $M$ inputs,
$\bm{X}=\bm{x}_{1}^{*},\ldots,\bm{x}_{M}^{*}$, sampled from the input space
$\mathcal{X}$. We generate a joint sample from the posterior at these
locations $(y_{1}^{*},\ldots,y_{M}^{*})\sim
p_{\bm{\theta}}(f(\bm{X})|\bm{y}_{1:n})$. We choose $\bm{x}_{n+1}$ to be the
input value with the lowest sampled value, $\bm{x}_{n+1}=\bm{x}_{j}^{\star}$
with $j=\arg\min_{i\in\\{1,...,M\\}}y_{i}^{*}$. If we want a batch of $q>1$
query points, we generate $q$ independent batches from the posterior over the
same set $\bm{X}$ and choose the lowest value from each batch (avoiding
duplicates).

###  3.2 Vecchia GP Approximations

Given $n$ GP observations $\bm{y}_{1:n}$ as described above, the Vecchia GP
approximation (Vecchia1988; Katzfuss2017a) uses a modified likelihood to turn
the $\mathcal{O}(n^{3})$ complexity for standard GP regression into
$\mathcal{O}(nm^{3})$, with $m\ll n$. (A further reduction to
$\mathcal{O}(nm^{2})$ is possible via more complicated grouping methods
described in Schafer2020.) The following provides details on Vecchia GPs,
including the likelihood and the posterior predictive distribution.

####  3.2.1 Vecchia Likelihood

Recall that any joint density of $n$ observations $\bm{y}_{1:n}$ can be
decomposed exactly as a product of conditional densities:
$p(\bm{y}_{1:n})=\prod_{i=1}^{n}p(y_{i}|\bm{y}_{1:i-1})$, where
$\bm{y}_{1:0}=\emptyset$. This motivates the Vecchia approximation
(Vecchia1988):

$\textstyle\widehat{p}(\bm{y}_{1:n})=\prod_{i=1}^{n}p(y_{i}|\bm{y}_{c(i)}),$
(2)

where each $c(i)\subset\\{1,\ldots,i-1\\}$ is a conditioning index set of size
at most $m$, and the $n$ conditional distributions $p(y_{i}|\bm{y}_{c(i)})$
are all Gaussian and can be computed in parallel, each in $\mathcal{O}(m^{3})$
time.

####  3.2.2 Vecchia Properties

Vecchia approximations have several attractive properties (e.g., as reviewed
by Katzfuss2020). The joint distribution
$\widehat{p}(\bm{y}_{1:n})=\mathcal{N}(\bm{\mu},\widehat{\bm{K}})$ implied by
(2) is multivariate Gaussian, where the inverse Cholesky factor
$\widehat{\bm{K}}^{-1/2}$ is sparse with $\mathcal{O}(nm)$ nonzero entries
(e.g., Katzfuss2017a). Under the sparsity pattern implied by the $c(i)$, the
Vecchia approximation results in the optimal $\widehat{\bm{K}}^{-1/2}$ as
measured by Kullback-Leibler (KL) divergence,
$KL(p(\bm{y})\|\widehat{p}(\bm{y}))$ (Schafer2020). Growing the conditioning
sets $c(i)$ decreases the KL divergence (Guinness2016a); for $m=n-1$, the
approximation becomes exact, $\widehat{p}(\bm{y})=p(\bm{y})$. Thus, $m$ trades
off low computational cost (small $m$) and high accuracy (large $m$);
crucially, high accuracy can be achieved even for $m\ll n$ in many settings.

####  3.2.3 Ordering and Conditioning Sets

Aside from the choice of $m$, the accuracy of a Vecchia approximation depends
on the ordering of $\bm{y}_{1:n}$ and on the choice of the conditioning sets
$c(i)$ in (2).

In our experience, the highest accuracy for a given $m$ can be achieved by
combining maximum-minimum-distance (maximin) ordering with conditioning on the
$m$ (previously ordered) nearest neighbors (NNs), which is illustrated in
Figure 1. Maximin sequentially picks each variable in the ordering as the one
that maximizes the minimum distance to previously ordered variables. Exact
maximin ordering and NN conditioning can be computed in quasilinear time in
$n$, and $m$ needs to grow only polylogarithmically with $n$ for
$\epsilon$-accurate Vecchia approximations with maximin ordering and NN
conditioning under certain regularity conditions (Schafer2020).

In Figure 1, we see that for index $i$ in the maximin ordering, the $m$
previously ordered NNs can be far away for small $i$ and very close for large
$i$. This global and local behavior helps Vecchia both learn accurate values
for hyperparameters and make good predictions. This is in contrast to low-rank
approximations that effectively use the same conditioning sets for each $i$
and hence focus on global behavior at the cost of ignoring local structure.

Figure 1: Illustration of maximin ordering and conditioning sets for $n=500$
uniformly sampled inputs (grey dots) on $[0,1]^{2}$. For $i=10$, $i=100$, and
$i=480$ in the left, middle, and right panel, we show the $i$th input (cyan
star) in the maximin ordering, the previously ordered $i-1$ inputs (blue
dots), including the (approximate) nearest $m=9$ inputs (red circles). For
increasing $i$, the first $i$ inputs form increasingly dense grids, so that
the distance to the $m$ NN decreases gradually and systematically. In
contrast, for fully independent conditional approximations (FIC), each input
has the same conditioning set, here illustrated by the first $m$ inputs in the
maximin ordering (magenta squares). Thus FIC effectively ignores many closer
previously-ordered inputs, which can lead to large approximation errors.

####  3.2.4 Vecchia Predictions

The Vecchia approximation can also be used to obtain an accurate approximation
of the joint posterior predictive distribution of an unobserved $n_{p}$-vector
$\bm{y}_{p}$ (Katzfuss2018) of the form

$\textstyle\widehat{p}(\bm{y}_{p}|\bm{y}_{1:n})=\prod_{i=1}^{n_{p}}p(y_{p}^{(i)}|\bm{y}_{c_{p}(i)})=\mathcal{N}(\bm{\mu}_{p},(\bm{L}^{T}\bm{L})^{-1}),$
(3)

where $\bm{y}_{c_{p}(i)}$ denotes the $m$ NNs of $y_{p}^{(i)}$ among
$\bm{y}_{1:n}$ (or $\bm{y}_{1:n}\cup\bm{y}_{p}^{(1:i-1)}$, if joint
predictions are desired). The distribution in (3) can be computed in
$\mathcal{O}(n_{p}m^{3})$ time, for example using the univariate Gaussian
conditional distributions in the middle part of (3).

Alternatively, the matrix expression for the posterior predictive distribution
on the right-hand side of (3) can be computed by first filling in the elements
of a sparse matrix $\bm{U}$, as follows. Let $g(i)$ be the indices of the
conditioning set for observation $i$ with input $\bm{x}_{i}$ and response
$y_{i}$. The $(j,i)$th element of $\bm{U}$ is given by:

$\bm{U}_{j,i}=\begin{cases}d_{i}^{-1/2},&i=j\\\ -b_{i}^{(j)}d_{i}^{-1/2},&j\in
g(i)\\\ 0,&\text{otherwise},\end{cases}$ (4)

where
$\textbf{b}_{i}^{T}=K(\bm{x}_{i},\bm{x}_{g(i)})K(\bm{x}_{g(i)},\bm{x}_{g(i)})^{-1}$,
$d_{i}=K(\bm{x}_{i},\bm{x}_{i})-\bm{b}_{i}^{T}K(\bm{x}_{g(i)},\bm{x}_{i})$. We
don’t need to construct the entire matrix $\bm{U}$, but instead we simply need
to fill in sub parts of the matrix for posterior inference. In particular let
$\bm{U}_{p,p}$ denote the sub-matrix of $\bm{U}$ corresponding to prediction
locations and $\bm{U}_{n,p}$ denote the sub-matrix for training locations and
prediction locations. We then use these sub-matrices to compute the posterior
mean,

$\displaystyle\bm{\mu}_{p}=(\bm{U}_{p,p}^{T})^{-1}(\bm{U}_{n,p})^{T}\bm{y}_{1:n},$

and the Cholesky factor of the posterior precision is simply given by
$\bm{L}=\bm{U}_{p,p}^{T}$. The matrices $\bm{U}_{p,p}$ and $\bm{U}_{n,p}$ are
both very sparse with at most $m$ non-zero elements per column, reducing the
cost of the operations. We provide additional details on Vecchia predictions
in Appendix A.2.1.

###  3.3 Vecchia and SVI-GP

Hensman2013GaussianData showed how inducing point GP approximations can be
trained using stochastic variational inference, and we refer to this model as
SVI-GP. To demonstrate how Vecchia GPs perform on functions with steep
regions, we compared Vecchia GP to SVI-GP on the Michalewicz function
(molga2005test).

We created training data by generating $n=500$ observations from the
15-dimensional version of the function and adding independent noise from
$\mathcal{N}(0,0.05^{2})$. We then generated a test set mimicking a typical BO
prediction setting by finding the minimum in the training set and generating
50 noiseless realizations in a neighborhood of the minimum.

For a range of values for $m$, we fitted an exact Gaussian process, the
Vecchia approximation and the SVI-GP, and computed all the posteriors on the
test set. Then we computed the KL-divergence
$KL(p(\bm{y}_{p}|\bm{y}_{1:n})\|\widehat{p}(\bm{y}_{p}|\bm{y}_{1:n}))$ between
the exact posterior and each of the approximate posteriors. We repeated this
process, including generating new data, and averaged the results as we varied
$m$. For Vecchia, $m$ corresponds to the number of NN, and for SVI-GP, $m$
corresponds to the number of inducing points.

Figure 2 shows the log of the KL-divergence for both models as we increase
$m$. In this example, Vecchia outperformed SVI-GP for all values of $m$ by
several orders of magnitude. Based on further exploration, we conjecture that
Vecchia more closely matches the local behavior of the function at the optima
and SVI-GP is more prone to return estimates similar to the average over the
entire domain. This success at modeling local fluctuations makes Vecchia a
strong candidate for BO, where we often have sharp spikes at the global
optima.

Figure 2: For prediction of the 15-dimensional Michalewicz function, the log
of the KL-divergence between an exact GP and approximate GPs based on the
Vecchia approximation (blue) and the SVI-GP approximation (red).

##  4 Vecchia GPs for Bayesian optimization

In this section we detail our modifications to Vecchia GPs to improve their
performance for BO.

###  4.1 Warped Kernel

A flexible GP kernel can be obtained by taking a base kernel $\tilde{K}$
(e.g., Matérn 5/2) and warping the inputs with a nonlinear function
$w(\cdot)$,

$K(\bm{x},\bm{x}^{\prime})=\tilde{K}(\|w(\bm{x})-w(\bm{x}^{\prime})\|).$ (5)

In particular, we warp using the Kumarswamy CDF (balandat2020botorch), which
reduces to the identity for a certain choice of hyperparameter values. To
fully utilize this flexibility within Vecchia GPs, we carry out the ordering
and conditioning-set selection based on Euclidean distance between the warped
inputs $w(\bm{x}_{1}),\ldots,w(\bm{x}_{n})$ instead of the original inputs
$\bm{x}_{1},\ldots,\bm{x}_{n}$. When the training data are observed with
noise, we assume $\tilde{K}$ includes a nugget term $\mathbbm{1}_{i=j}\tau$.

###  4.2 Training Using Stochastic Gradient Descent

The kernel $K=K_{\bm{\theta}}$ depends on unknown hyperparameters in the base
kernel $\tilde{K}$ and in the warping function $w$. To learn the potentially
high-dimensional hyperparameter vector $\bm{\theta}$, we optimize the Vecchia
loglikelihood $\log\widehat{p}_{\bm{\theta}}(\bm{y}_{1:n})$ given by (2) using
gradient descent. Importantly, on the log scale, this density is a sum of $n$
terms of the form $\log p(y_{i}|y_{c(i)})$, each of which can be computed in
$\order(m^{3})$ time (given the $m$ NNs). We use this to approximate the
gradient of (2) using mini-batches of the data:

$\textstyle\nabla_{\bm{\theta}}\widehat{p}(\bm{y}_{1:n})\approx\sum_{j\in
b}(n/|b|)\,\nabla_{\bm{\theta}}\log p(y_{j}|\bm{y}_{c(j)}),$ (6)

where $b\subset\\{1,...,n\\}$.

With this change, the training time depends on the mini-batch size $b$ instead
of the total data size $n$ given the ordering and NN. This reduction in time
complexity is critical, as we are estimating the GP hyperparameters after each
step in the BO procedure.

Figure 3 shows the negative log-likelihood vs the epoch (i.e., the number of
passes through the data) during training for different batch sizes, for data
generated using a simple sinusoid. This illustrates the typical training
dynamics for Vecchia GPs trained using mini-batch gradient descent. For of our
applications we find that a batch size of 64 offers a good balance between
stable training and speed.

Figure 3: Log-likelihood vs. epoch for different batch sizes.

###  4.3 Nearest Neighbors and Ordering

The nearest-neighbor search can become a bottleneck in the Vecchia procedure,
but we have found that using approximate NN yields good performance. In
particular, we suggest using approximate NNs with product quantization
(jegou2010product; JDH17). This allows for balancing accuracy and speed as we
consider how much data is available. As illustrated in Figure 1, in our
experience the approximate NN search is often good but not perfect for small
to medium maximin index $i$ (middle panel) but it is typically highly accurate
for large $i$.

For ordering data, we have found that an approximate maximin ordering works
well. The data is shuffled, split into disjoint sets, ordered within each set,
and recombined. Details of this procedure are given in Appendix 2. For really
large datasets, maximin ordering can even be replaced by a completely random
ordering at the cost of a slight decrease in accuracy for a given $m$
(Guinness2016a).

Together, these changes result in Vecchia approximations that can be applied
to 1,000-dimensional GPs with $n=100{,}000$ observations and beyond.

###  4.4 Calibration of the Prediction Variance

To improve uncertainty quantification for out-of-sample predictions, we
calibrate the marginal predictive variances as in Katzfuss2020. In short, we
maximize the log-probability of a hold-out set by adding a term to the
diagonal of the covariance matrix.

To do this we first obtain an estimate $\bm{\theta}^{*}$ of the
hyperparameters based on all available data. We then create a hold-out set,
$\tilde{\bm{y}}$, and compute the posterior mean
$\tilde{\bm{\mu}}_{\bm{\theta}^{*}}$ and variance
$\tilde{\bm{\Sigma}}_{\bm{\theta}^{*}}$ for this hold-out set based on the
training data not in the hold-out set, as described in Section 3.2.4. Finally,
we choose the variance inflation factor, $b_{v}$, to maximize the posterior
log-probability of this holdout set:

$\displaystyle\max_{b_{v}\in[0,2]}\text{ log
}p(\tilde{\bm{y}}|\tilde{\bm{\mu}}_{\bm{\theta}^{*}},\tilde{\bm{\Sigma}}_{\bm{\theta}^{*}}+b_{v}\bm{I}_{\tilde{n}}).$
(7)

The details of this procedure are given by Algorithm 1, where the funtion
$KNN(\bm{a},\bm{b},c)$ returns the $c$ NN of $\bm{b}$ in $\bm{c}$ and the
function $subsample(\bm{i},j)$ selects a subsample of size $j$ from $\bm{i}$.

With $b_{v}>0$, the predictive variance is inflated, which encourages more
exploration in the resulting BO procedure. Since the observations are
processed to have zero mean and variance equal to one, the value of $b_{v}$ is
generally at most 1.4 but usually less than 0.1. In our experience, this
variance correction is helpful when the observations are noisy or the target
function is highly variable. However, if the function is smooth and we have
noise-free observations, simply using the initial Vecchia predictions is
sufficient.

Algorithm 1 Vecchia Variance Calibration

Input: data $\textbf{x},\textbf{y}$, parameters $\bm{\theta}$

$I_{x^{*}}\leftarrow\arg\max(\textbf{y})$

$I_{O}\leftarrow KNN(\textbf{x}[I_{x^{*}}],\textbf{x},5*q)$

$I_{O}=subsample(I_{O},q)$

$\bm{\tilde{}}{X}=\textbf{x}[I_{O}]\cup\textbf{x}[I_{\bm{x}^{*}}]$

$\bm{\tilde{}}{y}=\textbf{x}[I_{O}]\cup\textbf{y}[I_{\bm{x}^{*}}]$

$\tilde{\mathcal{D}}=\\{\textbf{x}\setminus\bm{\tilde{}}{X},\textbf{y}\setminus\bm{\tilde{}}{y}\\}$

$\mathcal{L}_{b_{v}}=p(\bm{\tilde{}}{y};\tilde{\bm{\mu}}_{\bm{\theta}^{*}},\tilde{\bm{\Sigma}}_{\bm{\theta}^{*}})$

Return: $\arg\max_{b_{v}\in[0,2]}\mathcal{L}_{b_{v}}$

###  4.5 Speeding Up Existing BayesOpt Methods

Below we provide details for incorporating Vecchia GPs into existing BayesOpt
frameworks to improve their computational speed.

####  4.5.1 Vecchia with TuRBO

Trust Region Bayesian optimization (TuRBO) (eriksson2019scalable) is a
Bayesian optimization procedure that uses hyper-rectangles called trust
regions. In TuRBO, independent GPs are fit to data within the trust regions.
Query points are chosen using Thompson sampling, and the regions grow or
shrink based on how often they contain new optima. TuRBO allows for a single
trust region, TuRBO-1 or $t$ simultaneous trust regions, TuRBO-t. TuRBO-1 is
based on a single global GP fit to all data, whereas TuRBO-t fits $t$ local
GPs, each using only the data within one of the trust regions. The individual
trust regions help by allowing greater flexibility in modeling the local
behavior of the function and facilitating a more thorough exploration of the
entire input space (Eriksson2019).

Veccchia fits in the TuRBO framework without much change. Using either TuRBO-1
or TuRBO-t, the same (global) Vecchia GP is shared between all trust regions,
and predictions within each trust region are made using the $m$ NNs among all
data (i.e., the NNs are not restricted to be in the trust region for either
TuRBO-1 or TuRBO-t). By conditioning on the NNs, we emphasize the local
behavior of the target function, and by training on all the data we ensure our
hyperparameters are well calibrated. When using TuRBO-1, the variance
correction of Vecchia amounts to inflating or deflating the predictive
variance for points in and around the trust region. For the remainder of this
paper, we are referring to TuRBO-1 when we say TuRBO.

####  4.5.2 Vecchia With SAASBO

The Sparse AxisAligned Subspace BayesOpt (SAASBO) (eriksson2021high) places
horseshoe priors on the lengthscale parameters of the GP surrogate. Using a
standard exact GP, the time complexity to generate a sample from the posterior
is $\mathcal{O}(dn^{3})$. Vecchia GPs can improve on this by again reducing
the complexity required to invert the training data covariance matrix, as we
make the same NN approximation. We can also reduce the cost to evaluate the
likelihood further by using the data-splitting idea from neal2011mcmc. That is
we can compute the likelihood on batches of the data when we perform MCMC
sampling. Since BayesOpt queries points that can be close to one another, we
expect this procedure to work well, as this is the setting when neal2011mcmc
suggests data splitting ought to work well.

####  4.5.3 Vecchia With Parallel Acquisition Functions

For parallel acquisition functions such as q-EI or q-UCB, the query points
that maximize the acquisition function can be high-dimensional. Gradients of
these acquisition functions can ease this inner optimization problem, but the
closed-form expression may not be known. Fortunately, Monte Carlo estimates of
the gradient can be used (WilsonMaximizingOptimization). This approach allows
us to estimate the gradient of the utility function with respect to X (the
collection of query points) as:

$\nabla\ell(\textbf{y})=\frac{d\ell(\textbf{y})}{d\textbf{y}}\frac{d\textbf{y}}{d\mathcal{M}(\textbf{X)}}\frac{d\mathcal{M}(\textbf{X})}{d\textbf{X}},$
(8)

where $\mathcal{M}$ is our surrogate model. Under certain conditions, a Monte
Carlo approximation of (8) will be unbiased and will converge almost surely to
a set of stationary points of q-EI (wang2016parallel). In short, the following
can be used to estimate the gradient of q-EI in a stochastic optimization
algorithm and we should expect good performance:

$\nabla\mathcal{L}(\textbf{X})\approx\frac{1}{m}\sum_{k=1}^{m}\nabla\ell(\textbf{y}^{k}),$
(9)

where $y^{k}$ are draws from the posterior over $f$. We can then use (9) to
optimize parallel acquisition functions such as q-EI. Sufficient conditions
for the convergence of the stochastic optimization algorithm given in
wang2016parallel using the approximate gradient given by (9) requires, among
other things, that the mean and covariance functions of the GP be continuously
differentiable for each point in the domain. While this is not guaranteed to
be satisfied for Vecchia GPs at a set of points of measure zero, as the data
grows the effect of this boundary difference should diminish. In practice we
can optimize the acquisition function using standard gradient descent even for
small $n$. Additionally, since Vecchia prediction for joint samples is so fast
Vecchia GPs can greatly improve the speed of this inner optimization
especially for high dimensional problems.

##  5 Numerical Comparisons

We evaluated Vecchia GPs with TuRBO on a set of benchmarks. The synthetic
functions were Ackley-5, Hartmann-6 and Levy. For the Levy function we used
what we referred to as the Levy 55-20 function, where we warp the space with a
sigmoid type function and add irrelevent dimensions. Finally, we used two
reinforcement-learning examples found in eriksson2019scalable: Robot Pushing
(wang2018batched) and a Lunar Lander problem in OpenAI’s gym environment.

###  5.1 GP Surrogates

We compared the performance of the following GP surrogates within TuRBO-1:

Exact GP:

A standard Gaussian process, which is very slow for large $n$.

LR-First m:

Low-rank approximation where conditioning set taken as the first $m$ points in
the maximin ordering. The conditioning set is the same for all $n$ data
points.

SVI-GP:

The SVI GP (Hensman2013GaussianData) places a Gaussian variational
distribution on the outputs of “pseudo-points.” The variational parameters,
including the location of the $m$ pseudo points are learned with mini-batch
gradient descent.

Vecchia:

Our method as described in the previous sections, where the likelihood is
approximated by conditioning each observation on the $m$ NN subject to the
constraint that the neighbors appear before the data point in the maximin
ordering.

For the approximate methods, we let $m$ increase polylogarithmically with $n$,
as $m=C_{0}*\log^{2}_{10}(n)$. The value of $C_{0}$ was set to 7.2; for
details about this selection see Appendix A.3.2.

###  5.2 Comparison Measures

For target functions whose optimum is not known, we simply used the best value
found as our measure of performance. For functions whose optimum is known, we
used regret as our measure of performance. Regret after $t$ function
evaluations is given by, $r_{t}=\min_{1\leq i\leq t}\\{y_{i}-f^{*}\\}$, where
$f^{*}$ is the true minimum and the min is taken over all previously observed
function values. Our regret curves are averaged over independent repeats.

###  5.3 Results

####  5.3.1 Hartmann-6 and Ackley-5

Figure 4 shows the results for Hartmann-6 with batch size $q=20$. Vecchia with
variance inflation performed better than any other method, including the exact
GP. As this function is relatively smooth, in this case Vecchia benefits from
the ability to provide joint predictions and from variance inflation.

Figure 5 shows average log regret for Ackley-5 with a batch size of $q=20$.
For Vecchia, we did not use the variance correction. The results exceeded the
exact GP by 500 function evaluations and outperformed the other approximate
surrogates with ease.

Both these examples show the flexibility of the Vecchia GP and how modeling
choices can be made to improve performance on particular tasks.

Figure 4: Log regret for Hartmann-6 with batch size $q=20$. Curves are
averages of 30 repeats. Figure 5: Log regret for Ackley-5 with batch size
$q=20$. Curves are averages of 30 repeats.

####  5.3.2 Warped and embedded Levy

From Figure 6, we see that Vecchia (without warping or variance inflation) did
better than all other surrogates; by including the warping capability, the
performance of Vecchia improved further. We believe that the warping did not
improve the exact GP the same way it improved Vecchia because Vecchia with
warping is making the best local approximations of the function possible. Due
to the NN conditioning, Vecchia was able to focus on local behavior of this
relatively complex function and ignored irrelevant distant parts of the
function.

Figure 6: Log regret for the warped and embedded 20-55 Levy function with
batch size $q=50$. Curves were averaged over 30 repeats.

####  5.3.3 Robot Pushing

In the $14$-dimensional robot-pushing problem (wang2018batched), there are two
robot arms pushing objects towards a target location. We must minimize the
distance to the target location.

As shown in Figure 7, Vecchia performed the best among all methods, eventually
even outperforming the exact GP for a large number of function evaluations. In
this setting, the response was noisy and the target function is highly
variable, which is a situation in which using variance inflation for Vecchia
GPs within BO works best.

Figure 7: Best value found for 14D Robot Pushing problem with batch size
$q=50$. Curves were averaged over 20 repeats.

####  5.3.4 Lunar Lander

We must tune a PID controller (eriksson2019scalable) to successfully land a
Lunar module at a fixed point on a set of 50 maps. The results are noisy,
because each run of the function is an average over 50 experiments.

Figure 8 shows that scaled Vecchia and SVI-GP performed best, both
outperforming the exact GP. Just as with the 14D robot pushing problem, the
variance inflation was helpful because the function is noisy.

Figure 8: Best value found for 12$D$ Lunar Lander problem with batch size
$q=50$. Curves are averaged over 20 repeats.

##  6 Discussion

We explored the use of approximate Vecchia GPs for speeding up Bayesian
optimization (BO). We described how training of Vecchia GPs can be made
efficient through mini-batch gradient descent, approximate nearest neighbors
and approximate variable orderings. Additionally, we incorporated Vecchia GPs
within the TuRBO BO procedure, and commented on how to use Vecchia GPs for
other BO frameworks such as SAASBO and parallel acquisition functions.

In our numerical experiments on a suite of test functions, we found Vecchia
GPs to often be as or even more accurate than other GP surrogates, including
the exact GP. This is remarkable, in that, in addition, our Vecchia
approximations drastically reduce the computational complexity relative to the
exact GP, with a computational cost that largely depends (linearly) on the
mini-batch size. As a result, Vecchia can be even computationally cheaper than
approximate pseudo-point methods such as SVI-GP, which must optimize over a
large number of hyperparameters related to the pseudo-points, especially in
high-dimensional input spaces, while Vecchia does not reduce any additional
hyperparameters beyond those already present in the exact-GP model to be
approximated. We conjecture that Vecchia’s accuracy improvements are related
to the variance correction and the global-local nature of Vecchia GPs; the
nearest neighbor portion of Vecchia focuses on the local behavior of the
target function during prediction while using the entire dataset to estimate
the hyperparameters from the global information. Together, Vecchia GPs
naturally balances local and global information when estimating the target
function.

Our results suggest that Vecchia GPs can be useful for speeding up BO when
many evaluations of the target function are feasible and necessary, and we are
hopeful that Vecchia approximations can be a standard tool for BO in such
settings.

There are several potential avenues for future work. Vecchia GPs are
competitive without GPU computing and thus enable fast BO even without access
to sophisticated computational resources; however, Vecchia GPs are in
principle well suited for GPU-type computations and can thus be sped up
further. For highly noisy target functions, we believe that further accuracy
gains can be made by latent Vecchia approximations with incomplete Cholesky
decomposition (Schafer2020) within Thompson sampling. Further, correlation-
based Vecchia GPs (Kang2021) may be well suited for neural architecture search
(kandasamy2018neural).

##  Appendix A Appendix

###  A.1 Approximate Variable Orderings

The function below gives an approximate maxmin ordering which works well
empirically, and can be run in parallel. We assume the inputs have been
processed to the correct space (e.g. warped, scaled by inverse lengthscales).
Additionally, we assume there is a function that can perform the exact maxmin
ordering, such as Algorithm C.1 in the supplementary material of
(Schafer2020).

Algorithm 2 Approximate MaxMin Ordering

1: Input: data x, parameters $MaxSubsetSize$

2: $N\leftarrow length(\textbf{x})$

3: if $N>MaxSubsetSize$ then

4: $order1\leftarrow
ApproximateMaxMin(\textbf{x}[0:floor(N/2)],MaxSubsetSize)$ {Run in Parallel.}

5: $order2\leftarrow
ApproximateMaxMin(\textbf{x}[floor(N/2):N],MaxSubsetSize)$ {Run in Parallel.}

6: $order=concatenate(order1,order2)$

7: else

8: $order=GetExactMaxMin(\bm{x})$

9: end ifReturn $order$

###  A.2 Prediction

Following Katzfuss2018 we consider the response first ordering (RF) using RF-
standard (RF-stand). In RF-stand, the $j$th prediction location conditions on
any previously ordered data point, including another prediction location. This
means we take into account the posterior covariance between some prediction
locations if they are closer together than to training data. For clarity we
use the subscribt $r$ to refer to the observed data and the subcript $p$ to
refer to predictions.

For prediction we form two matrices $\bm{U}_{r,p}$ and $\bm{U}_{p,p}$ whose
entries can be computed using equation 12. Using the $\bm{U}$ matrices we can
get the posterior mean and the Cholesky factor of the posterior variance with
complexity dependent on the number of prediction locations, and not the sample
size of the training data. See Section 4.1 in Katzfuss2018 for details.

The Cholesky factor of the posterior precision for prediction locations is
given by:

$\bm{L}_{p}=(\bm{U}_{p,p}^{-1})^{T}$ (10)

The posterior predictive mean is given by

$\bm{\mu}_{p}=(\bm{U}_{p,p}^{T})^{-1}\bm{U}_{r,p}^{T}\bm{y}_{1:n}.$ (11)

####  A.2.1 U Matrix

We before consider the matrix U, whose elements are given by:

$\textbf{U}_{j,i}=\begin{cases}d_{i}^{-1/2},&i=j\\\
-b_{i}^{(j)}d_{i}^{-1/2},&j\in g(i)\\\ 0,&\text{otherwise},\end{cases}$ (12)

When we allow the prediction locations to condition on other prediction
locations then there will be non-zero elements in the matrix
$\textbf{U}_{pp}$, and there will be columns with less than $m$ non-zero
elements in the matrix $\textbf{U}_{n,p}$. The significance of this becomes
obvious when we consider the posterior mean, since a location may condition on
less than $m$ observations. The prediction will instead rely on the value that
was predicted at other unknown locations.

This works well when the number of observations is much greater than the
number of predictions. If we restrict ourselves to condition only on observed
data then we ignore the joint nature of the posterior, but gain speed as all
the operations can be done in parallel.

####  A.2.2 Kernel and Warping

When the target function cannot be modeled by a stationary GP, it can help to
first warp the input space through a non-linear invertible function. To
incorporate this tool into the Vecchia framework, we consider a general kernel
that allows for warping,

$K(\bm{x},\bm{x}^{\prime})=\tilde{K}(\|w(\bm{x})-w(\bm{x}^{\prime})\|),$ (13)

x where $\|\cdot\|$ denotes Euclidean distance and $\tilde{K}$ is an isotropic
covariance function. We assume the warping function to operate element-wise on
each of the $d$ input dimensions, such that
$w(\bm{x})=\big{(}w_{1}(\bm{x}_{1}),\ldots,w_{d}(\bm{x}_{d}))^{\top}$, where
$w_{j}(\bm{x}_{j})=\omega_{\bm{\phi}_{j}}(\bm{x}_{j})/\lambda_{j}$; here,
$\lambda_{j}$ is the range parameter for the $j$th input dimension, and
$\omega_{\bm{\phi}}$ is bijective, continuous, and differentiable with respect
to $\bm{\phi}$. For our numerical experiments, we used the Kumarswamy CDF
(balandat2020botorch), $\omega_{a,b}(\bm{x})=1-(1-\bm{x}^{a})^{b}$; when its
parameters are set to $a=b=1$, (13) becomes the well-known automatic relevance
determination (ARD) kernel.

For approximating a GP with a warping kernel, we propose the warped Vecchia
approximation. Specifically, in $\eqref{eq:vecchia}$, we carry out the maximin
ordering of $\bm{y}_{1:n}$ and the selection of NN conditioning sets $c(i)$
based on Euclidean distance between the corresponding warped inputs,
$w(\bm{x}_{1}),\ldots,w(\bm{x}_{n})$. For prediction of a new variable
$y^{*}_{i}$ at $\bm{x}_{i}^{*}$, we also select $c^{*}(i)$ as the $m$ NNs to
$w(\bm{x}_{i}^{*})$ among $w(\bm{x}_{1}),\ldots,w(\bm{x}_{n})$. In the special
case of linear warping with $\omega(\bm{x})=\bm{x}$, this warped Vecchia
approach is equivalent to the scaled Vecchia approximation for emulating
computer experiments (Katzfuss2020).

###  A.3 Experiment Details

Below are extra details on the experimental setup.

####  A.3.1 Experimental Setup

During every round, each GP surrogate model under consideration was fit using
all of the available data. Thompson sampling was used within each TuRBO trust
region, with the input locations chosen from a Sobol sequence. We evaluated
the GP at $\min(5000,\max(5000,100d))$ locations $q$-separate times. We chose
the batch of $q$ elements to be the locations with the $q$ highest simulated
values within Thompson sampling. For all surrogates, we used a Matérn2.5
covariance function for $\tilde{K}$ in (13), assuming linear warping (i.e.,
$a=b=1$) for all examples except the Levy function. For approximate surrogate
models, we chose $m$ as described in Section A.3.2. The input space for all
test functions and problems was mapped to $[0,1]^{d}$. All results were
averaged over independent replicates, with starting points again chosen using
a Sobol sequence.

####  A.3.2 Choosing $m$

For the approximate surrogates, we grow $m$ polylogarithmically with $n$,
$m=C_{0}*\log^{2}_{10}(n)$. To simplify our comparison between surrogate
models, we chose a single value of $C_{0}$ for all our experiments. To choose
this $m$, we looked at the regret vs number of samples when we used the
Vecchia approximation with EI as our acquisition function and compared with
the regret curve using the exact GP. When using the Six-Hump camel function we
found a $C_{0}$ value of around 7.2 is the lowest for which Vecchia tracked
the exact solution closely. We excluded the Six-Hump camel function from the
comparison of regret curves to put the competing approximations on more even
footing.

####  A.3.3 Levy 55-20 function

We used a modified version of the Levy-20 function to show the benefit of
warping with Vecchia. We call this function Levy 55-20, and it is simply the
Levy-20 function with 35 extra irrelevant dimensions. The hyper cube on
$[0,1]^{55}$ was warped using $w(\bm{x})=S(4*\bm{x}-1)$ ($S(\cdot)$ is the
sigmoid function) before being scaled to $[-5,5]^{55}$. Only the first 20
dimensions are used to evaluate the Levy-20 function.