††thanks: AID funding

# Experimental investigation of the thermal emission cross-section of nano-
resonators using hierarchical Poisson-disk distributions

Denis Langevin  DOTA, ONERA, Université Paris-Saclay, F-91123 Palaiseau -
France  Clément Verlhac  DOTA, ONERA, Université Paris-Saclay, F-91123
Palaiseau - France  Julien Jaeck  DOTA, ONERA, Université Paris-Saclay,
F-91123 Palaiseau - France  Loubnan Abou-Hamdan  DOTA, ONERA, Université
Paris-Saclay, F-91123 Palaiseau - France  Institut Langevin, ESPCI Paris, PSL
University, CNRS, 1 rue Jussieu, F-75005 Paris, France  Nathalie Bardou
Center for Nanosciences and Nanotechnology (C2N) - CNRS, Université Paris-
Saclay, 10 Boulevard Thomas Gobert, 91120 Palaiseau, France  Christophe Dupuis
Center for Nanosciences and Nanotechnology (C2N) - CNRS, Université Paris-
Saclay, 10 Boulevard Thomas Gobert, 91120 Palaiseau, France  Yannick De Wilde
Institut Langevin, ESPCI Paris, PSL University, CNRS, 1 rue Jussieu, F-75005
Paris, France  Riad Haïdar  DOTA, ONERA, Université Paris-Saclay, F-91123
Palaiseau - France  Patrick Bouchon<EMAIL_ADDRESS>DOTA, ONERA,
Université Paris-Saclay, F-91123 Palaiseau - France

###### Abstract

Effective cross-sections of nano-objects are fundamental properties that
determine their ability to interact with light. However, measuring them for
individual resonators directly and quantitatively remains challenging,
particularly because of the very low signals involved. Here, we experimentally
measure the thermal emission cross-section of metal-insulator-metal nano-
resonators using a stealthy hyperuniform distribution based on a hierarchical
Poisson-disk algorithm. In such distributions, there are no long-range
interactions between antennas, and we show that the light emitted by the
metasurface behaves as the sum of cross-sections of independent nanoantennas,
enabling direct retrieval of the single resonator contribution. The emission
cross-section at resonance is found to be of the order of
$\mathbf{\lambda_{0}^{2}/3}$, a value that is nearly three times larger than
the theroretical maximal absorption cross-section of a single particle but
remains smaller than the maximal extinction cross-section. This measurement
technique can be generalized to any single resonator cross-section, and we
also apply it here to the extinction cross-section.

Optical nanoantennas, often inspired by radiofrequency designs of classical
antennas [1], have been used to tame electromagnetic wave properties at
optical wavelengths far below the diffraction limit [2, 3, 4, 5, 6]. The
interactions of single nanoantennas with light can be of great interest to a
multitude of light applications [7, 8, 9, 10, 11, 12, 13, 14], but they can
also be used as building blocks (meta-atoms) for larger metasurfaces [15, 16].
Among the fundamental properties of single objects are their effective
absorption and scattering cross-sections, the sum of which is the extinction
cross-section [17, 18]. These quantities describe the single object’s ability
to interact with light at a given frequency but are often complex to quantify
experimentally [19, 20].

In practice, the direct measurement of the absorption cross-section has so far
remained rather elusive [21]. As a result, the absorption cross-section is
typically obtained by taking the difference between the extinction and
scattering cross-sections. The emission cross-section can also be introduced
based on the local polarized emissivity density that is linked to absorptivity
by Kirchhoff’s law [22]. Direct measurement of infrared thermal emission of
single objects has been previously demonstrated but the radiated signal is
very low, particularly in comparison with the emitting background [23, 24].
Determining absorption and emission cross-sections of single objects is of the
utmost importance for applications, such as thermoplasmonics [25], energy
conversion, [26, 27] and meta-emitters with tailored emission (polarization,
spectral) [28, 29, 30, 31, 32].

In this letter, we measure, both directly and quantitatively, the thermal
emission cross-section of nanoresonators arranged as a correlated disorder
distribution. The metasurface is built following a hierarchical Poisson disk
algorithm with various densities that gives a stealthy hyperuniform
distribution. For low densities, we show that the thermal emission is linearly
dependent on the density of nanoantennas, thus behaving as an independent
combination. In that case, we measure the sum of the individual antennas’
thermal emission cross-sections. These measurements are shown to be in fair
agreement with electromagnetic computations on a single nanoantenna.
Eventually, we demonstrate that our experimental protocol is easily adaptable
to other cross-section measurements by applying it to the measurement of
extinction cross-sections.

Figure 1:  Illustration of the (a) scattering, (b) absorption and (c) emission
cross-sections for a nano-object. (d) Non-overlapping cross-sections (e)
Overlapping cross-sections.

The various electromagnetic cross-sections of a single nano-object on an
opaque substrate are illustrated in Fig. 1(a-c). The scattering cross-section
$\sigma_{scat}$ and the absorption cross-section $\sigma_{abs}$ are shown for
the case of a normally incident monochromatic plane wave. Their sum gives the
well-known extinction cross-section $\sigma_{ext}=\sigma_{abs}+\sigma_{scat}$.
At a temperature T, the nano-object emits light due to thermal fluctuations
according to Planck’s law. Thus, we introduce the emission cross-section
$\sigma_{em}$ that describes its ability to emit light for a given wavelength
$\lambda$, direction of emission $\vec{u}$ and polarization. The local form of
Kirchhof’s law equates the absorption cross-section to the emission cross-
section for similar conditions, i.e.,
$\sigma_{em}(\vec{u},\lambda)=\sigma_{abs}(\vec{u},\lambda)$[22].

An object’s emission cross-section in a given direction $\vec{u}$ at a
wavelength $\lambda$ is defined by:

$\sigma_{em}(\vec{u},\lambda)=\int_{V}\eta(\vec{u},\vec{r},\lambda)d^{3}\vec{r}$
(1)

where $\eta(\vec{u},\vec{r},\lambda)$ is the local spectral directional
emissivity at position $\vec{r}$ in the object. This is equal to the
absorptivity of the object and can be computed using the object’s local
permitivitty [22].

Figure 1(d-e) extends this illustration to assemblies of nano-objects emitting
light. In the first case, the nano-objects are sufficiently far apart so that
no interactions may arise. The total emission cross-section
$\sigma_{em}^{(tot)}$ is then given by the sum of the individual emission
cross-sections $\sum\sigma_{em}^{(res)}$. In the second case (Fig. 1(e)), the
resonators are much closer to one another and can no longer be considered
independent. Therefore, the individual cross-sections become spatially
overlapping, giving rise to near-field interactions so that
$\sigma_{em}^{(tot)}\neq\sum\sigma_{em}^{(res)}$. In most situations, the
overlap of cross-sections leads to a loss of efficiency, i.e.,
$\sigma_{em}^{(tot)}<\sum\sigma_{em}^{(res)}$. However, hybridized modes may
appear and generate constructive interferences, leading to
$\sigma_{em}^{(tot)}>\sum\sigma_{em}^{(res)}$. In any case, the total far-
field emission obviously remains below the fundamental limit set by the fact
that the emissivity of any object cannot be larger than 1.

Periodic arrays of nanoresonators are commonly employed to reach this
theoretical limit, in which complete absorption of incoming light is readily
achieved by reaching a critical coupling condition [33].

However, the array period, along with diffraction effects, play a major role
in the overall behavior of the sample, concealing the contribution of the
individual resonators [34]. The emission cross-section of a nanoantenna in an
array of period $d$ at critical coupling can be written as
$\sigma_{em}^{(tot)}=d^{2}$ at resonance. Consequently, the measurement of
electromagnetic cross-sections can only be effected in the case of a single
resonator or an aperiodic resonator distribution.

In what follows, the emission response of hyperuniform disordered metasurfaces
using metal-insulator-metal (MIM) resonators as building blocks is
investigated, as depicted in Fig. 2. This type of antenna, commonly used in
metasurface design, acts as a Fabry-Perot nano-cavity [33, 35]. Here, the MIM
resonators consist of square gold antennas placed on top of a zinc sulfide
(ZnS) layer with a metallic backplate. The thickness of the dielectric layer
is 280 nm, and the antenna width $w$ is chosen so that the resonance
wavelength is in the long-wave infrared spectral range (8-12
$\mathrm{\SIUnitSymbolMicro m}$). The metallic layer thickness is chosen to be
significantly larger than the skin-depth of the metal. A scheme of this
resonator structure is shown in Fig. 2(a-c), together with an SEM image of a
single resonator.

Figure 2: (a) Scheme of the sample, with the overall distribution. (b) Scheme
of the single resonator geometry. (c) Scanning electron microscope (SEM) top
view image of a single resonator ($w=1.6\,$\mathrm{\SIUnitSymbolMicro m}$$).
(d) Structure factor of a distribution where $N=40000$ and $L=1$ mm. (e)
Examples of three consecutive states of the resonator distribution, after N
steps. The dashed circles are exclusion zones where no additional antennas can
be placed.

The stealthy hyperuniform disordered distribution is obtained following a
hierarchical Poisson-disk algorithm [36]. This algorithm resembles Random
Sequential Adsorption (RSA) [37], in that it builds the distribution gradually
by randomly picking new positions one by one, only adding a new position to
the distribution if it satisfies a pre-defined separation distance with
respect to all previously chosen positions. The hierarchical Poisson-disk
algorithm introduces a large exclusion distance between new positions, that is
gradually decreased as the distribution density $D$ increases.

This ensures that the algorithm can place a precisely chosen number of
resonators $N$, rather than stopping when the exclusion zones overlap, and
that the $n$ chosen positions first are also hyperuniform, for any $n<N$.
Thus, once a design of given density $D$ is created, designs of any density
$d<D$ can be retrieved by simply taking a subset of the densest distribution.
This guarantees that the density variations between samples are only caused by
adding or removing resonators from an existing distribution, without changing
the positions of other resonators in the process. This type of distribution
minimizes the long-range correlations between resonators, thereby precluding
the contribution of structural effects of the distribution on the signal [38,
39]. A way to characterize hyperuniformity is through the structure factor of
the resonator distribution that can be computed for any distribution
$\\{R\\}_{j}$ of $N$ points as [39]:

$S(\vec{k})=\frac{1}{N}\sum_{j_{1}}\sum_{j_{2}}\exp^{-i\vec{k}(\vec{R_{j_{1}}}-\vec{R_{j_{2}}})}.$

The structure factor for a hierarchical Poisson-disk distribution with
$N=40000$ is plotted in Fig. 2(d). The dark central area in this structure
factor is a typical characteristic of a stealthy hyperuniform distribution
[40]. The absence of high intensity peaks outside of the dark area indicates
that there are no long-range correlations between positions. An example of the
state of the distribution given by the hierarchical Poisson-disk algorithm at
different steps is shown in Fig. 2(e). The overall resonator density $D$ of
each sample is determined by the number of resonators per sample area:
$D=N/L^{2}$. We can identify samples by their resonator density $D$, or by the
period $L_{eq}=1/\sqrt{D}$ of its equivalent periodic array of the same
density.

Figure 3: (a) Measured emissivity $\epsilon$ spectra of the hyperuniform
disordered metasurface with resonator density
$D=1/64\,$\mathrm{\SIUnitSymbolMicro m}$^{-2}$, or equivalent period
$L_{eq}=8\,$\mathrm{\SIUnitSymbolMicro m}$$ (blue curve) and
$D=1/256\,$\mathrm{\SIUnitSymbolMicro m}$^{-2}$, or
$L_{eq}=16\,$\mathrm{\SIUnitSymbolMicro m}$$ (red curve). The spectra were
measured at $T=523\,\mathrm{K}$. The SEM images of the two metasurfaces are
shown above at the same magnification. (b) Emissivity as a function of the
density of the metasurface at resonance wavelength
$\lambda_{r}=8.7\,$\mathrm{\SIUnitSymbolMicro m}$$ for several measurements
(blue dots). The linear regression slope (orange line) gives the emission
cross-section $\sigma_{em}$ at $\lambda_{r}=8.7\,$\mathrm{\SIUnitSymbolMicro
m}$$.

The infrared thermal emission of disordered MIM metasurfaces with various
densities, heated at $T=523\,\mathrm{K}$ with a Linkam THMS600 heating stage,
is collected around normal incidence by a parabolic mirror. The heated
metasurface acts as an external source for a Fourier transform infrared
spectrometer (FTIR) that measures its emission spectrum.

The studied metasurfaces consist of square-patch MIM antennas of width
$w=1.6\,$\mathrm{\SIUnitSymbolMicro m}$$ distributed according to the
hierarchical Poisson-disk distribution at various densities. The emissivity
spectra of metasurfaces with resonator density
$D=1/64\,$\mathrm{\SIUnitSymbolMicro m}$^{-2}$ and
$D=1/256\,$\mathrm{\SIUnitSymbolMicro m}$^{-2}$ are shown in Fig. 3(a), with
scanning electron microscope (SEM) images illustrating the distribution. Both
metasurfaces exhibit a resonance at the same wavelength
$\lambda_{r}=8.7\,$\mathrm{\SIUnitSymbolMicro m}$$, and as expected, the
measured emissivity for the denser metasurface is higher.

At each wavelength $\lambda$, the emission cross-section
$\sigma_{em}(\lambda)$ can be extracted from the slope of the linear
regression of the emissivity as a function of density. Indeed, as illustrated
in Fig. 1(d), for independant nano-emitters, the emissivity is given by
$\epsilon=\sigma_{em}D$. This behavior is illustrated at the resonance
wavelength in Fig. 3(b). The measured value gives
$\sigma_{em}=24\,$\mathrm{\SIUnitSymbolMicro m}$^{2}$. This means that
$\sigma_{em}\sim\frac{\lambda_{r}^{2}}{3}$ for the resonant wavelength of 8.7
$\mathrm{\SIUnitSymbolMicro m}$. This order of magnitude is to be expected for
a nanoresonator on a reflective backplate [19, 41], but this value is higher
than the maximum theoretical value for a single resonant dipolar particle,
i.e., $\sigma_{\mathrm{abs}}^{\mathrm{dipole}}=\frac{3\lambda^{2}}{8\pi}$ [42,
43].

Figure 4: Comparison of measured and simulated emission cross-section spectra.
The measurements correspond to a resonator distributions of density
$L_{eq}=16\,$\mathrm{\SIUnitSymbolMicro m}$$, while the electromagnetic
simulations are for single resonators using perfectly matched layers.

Following the same process at each wavelength, we obtain the emission cross-
section spectra, which are are shown on Fig. 4 for various resonator
geometries. This shows that both the resonance wavelength and the maximum
cross-section increase with the resonator width $w$, as expected [24].

To further validate these experimental measurements, we performed finite
element method (FEM) electromagnetic simulations (COMSOL Multiphysics [44]) of
single resonators using perfectly matched layers. The emission cross-section
of an individual square MIM resonator was retrieved and compared to the
experimental measurements in Fig. 4. The FEM simulations show very good
agreement with the experimental measurements, both in resonance wavelength and
magnitude.

Figure 5: (a) Measured reflectivity spectra for a resonator density
$D=1/25\,$\mathrm{\SIUnitSymbolMicro m}$^{-2}$, or equivalent period
$L_{eq}=5\,$\mathrm{\SIUnitSymbolMicro m}$$ (blue),
$D=1/64\,$\mathrm{\SIUnitSymbolMicro m}$^{-2}$, or
$L_{eq}=8\,$\mathrm{\SIUnitSymbolMicro m}$$ (red) and
$D=1/144\,$\mathrm{\SIUnitSymbolMicro m}$^{-2}$, or
$L_{eq}=12\,$\mathrm{\SIUnitSymbolMicro m}$$ (green). The measurements are
shown for an incidence angle $\theta=13^{\circ}$ (see inset). (b) Linear
regression of measured extinction against density to compute the extinction
cross-section.

We underscore here that the measurement protocol that we have presented, using
stealthy hyperuniform samples of varying antenna densities, is extremely
robust and can be generalized to any kind of cross-section. Concordantly, we
have also used it to measure extinction cross-sections. In this case,
reflectivity measurements were performed at an angle of incidence
$\theta=13^{\circ}$. Figure 5(a) shows the measured reflectivity spectra for
three different sample densities, showing a dip in reflectivity at the
resonance, representing the extinction by the individual nanoresonators that
increases with increasing density. The extinction values measured as a
function of density (Fig. 5(b)) follow a clear linear behavior at low
densities. This leads to $\sigma_{ext}=30\,$\mathrm{\SIUnitSymbolMicro
m}$^{2}$, a value close to $\frac{2\lambda^{2}}{5}$. Contrary to the emission
cross-section, this value is lower than the maximum theoretical value for a
dipolar particle
($\sigma_{\mathrm{ext}}^{\mathrm{dipole}}=\frac{3\lambda^{2}}{2\pi}$).
Saturation can be observed at higher densities, where the measured extinction
is no longer proportional to the density of resonators. This behavior
originates from the overlap between resonator cross-sections, decreasing the
efficiency of each individual resonator. In the current study, this was only
observed for the extinction measurements, as the resulting cross-section is
larger than that obtained from emissivity measurements. The highest sample
densities we characterized were sufficient to result in an extinction cross-
section overlap but were still too low to give rise to an emission cross-
section overlap.

In summary, we have presented a protocol to quantitatively measure single
resonator cross-sections using large distribution far-field measurements. We
applied this protocol to emission cross-sections of MIM resonator
distributions, for which we obtained
$\sigma_{em}=24\,$\mathrm{\SIUnitSymbolMicro m}$^{2}$. These results were
found to be in fair agreement with FEM simulations of single resonators,
confirming that our measurement protocol gives direct access to the single
resonator emission cross-section. Finally, we extended our measurement
protocol to extinction cross-sections. In this case, a saturation behavior was
observed at high densities, which stresses the importance of performing
measurements at various sample densities in order to properly isolate the
linear regime, in which the presented protocol is applicable.

## References

  * Munk [2003] B. A. Munk,  _Finite antenna arrays and FSS_ (John Wiley & Sons, 2003).
  * Brongersma [2008] M. L. Brongersma, Engineering optical nanoantennas, Nature Photonics 2, 270 (2008).
  * Novotny and Van Hulst [2011] L. Novotny and N. Van Hulst, Antennas for light, Nature Photonics 5, 83 (2011).
  * Biagioni  _et al._ [2012] P. Biagioni, J.-S. Huang, and B. Hecht, Nanoantennas for visible and infrared radiation, Reports on Progress in Physics 75, 024402 (2012).
  * Krasnok  _et al._ [2013] A. E. Krasnok, I. S. Maksymov, A. I. Denisyuk, P. A. Belov, A. E. Miroshnichenko, C. R. Simovski, and Y. S. Kivshar, Optical nanoantennas, Physics-Uspekhi 56, 539 (2013).
  * Koenderink [2017] A. F. Koenderink, Single-photon nanoantennas, ACS Photonics 4, 710 (2017).
  * Neubrech  _et al._ [2008] F. Neubrech, A. Pucci, T. W. Cornelius, S. Karim, A. García-Etxarri, and J. Aizpurua, Resonant plasmonic and vibrational coupling in a tailored nanoantenna for infrared detection, Physical Review Letters 101, 157403 (2008).
  * Brongersma  _et al._ [2015] M. L. Brongersma, N. J. Halas, and P. Nordlander, Plasmon-induced hot carrier science and technology, Nature Nanotechnology 10, 25 (2015).
  * Mann  _et al._ [2016] S. A. Mann, S. Z. Oener, A. Cavalli, J. E. Haverkort, E. P. Bakkers, and E. C. Garnett, Quantifying losses and thermodynamic limits in nanophotonic solar cells, Nature Nanotechnology 11, 1071 (2016).
  * Taylor and Zijlstra [2017] A. B. Taylor and P. Zijlstra, Single-molecule plasmon sensing: current status and future prospects, ACS Sensors 2, 1103 (2017).
  * Zhu  _et al._ [2019] L. Zhu, A. Fiorino, D. Thompson, R. Mittapally, E. Meyhofer, and P. Reddy, Near-field photonic cooling through control of the chemical potential of photons, Nature 566, 239 (2019).
  * Gillibert  _et al._ [2020] R. Gillibert, M. Malerba, D. Spirito, V. Giliberti, L. Li, A. G. Davies, E. H. Linfield, L. Baldassarre, R. Colombelli, and M. Ortolani, Nanospectroscopy of a single patch antenna strongly coupled to a mid-infrared intersubband transition in a quantum well, Applied Physics Letters 117, 101104 (2020).
  * Sakat  _et al._ [2020] E. Sakat, L. Wojszvzyk, J.-J. Greffet, J.-P. Hugonin, and C. Sauvan, Enhancing light absorption in a nanovolume with a nanoantenna: theory and figure of merit, ACS Photonics 7, 1523 (2020).
  * Rajabali  _et al._ [2022] S. Rajabali, S. Markmann, E. Jöchl, M. Beck, C. A. Lehner, W. Wegscheider, J. Faist, and G. Scalari, An ultrastrongly coupled single terahertz meta-atom, Nature Communications 13, 1 (2022).
  * Meinzer  _et al._ [2014] N. Meinzer, W. L. Barnes, and I. R. Hooper, Plasmonic meta-atoms and metasurfaces, Nature Photonics 8, 889 (2014).
  * Lalanne and Chavel [2017] P. Lalanne and P. Chavel, Metalenses at visible wavelengths: past, present, perspectives, Laser & Photonics Reviews 11, 1600295 (2017).
  * Bohren and Huffman [2008] C. F. Bohren and D. R. Huffman,  _Absorption and scattering of light by small particles_ (John Wiley & Sons, 2008).
  * Staude  _et al._ [2013] I. Staude, A. E. Miroshnichenko, M. Decker, N. T. Fofang, S. Liu, E. Gonzales, J. Dominguez, T. S. Luk, D. N. Neshev, I. Brener, _et al._ , Tailoring directional scattering through magnetic and electric resonances in subwavelength silicon nanodisks, ACS Nano 7, 7824 (2013).
  * Husnik  _et al._ [2012] M. Husnik, S. Linden, R. Diehl, J. Niegemann, K. Busch, and M. Wegener, Quantitative experimental determination of scattering and absorption cross-section spectra of individual optical metallic nanoantennas, Physical Review Letters 109, 233902 (2012).
  * Crut  _et al._ [2014] A. Crut, P. Maioli, N. Del Fatti, and F. Vallée, Optical absorption and scattering spectroscopies of single nano-objects, Chemical Society Reviews 43, 3921 (2014).
  * Baffou  _et al._ [2012] G. Baffou, P. Bon, J. Savatier, J. Polleux, M. Zhu, M. Merlin, H. Rigneault, and S. Monneret, Thermal imaging of nanostructures by quantitative optical phase analysis, ACS Nano 6, 2452 (2012).
  * Greffet  _et al._ [2018] J.-J. Greffet, P. Bouchon, G. Brucoli, and F. Marquier, Light emission by nonequilibrium bodies: local kirchhoff law, Physical Review X 8, 021008 (2018).
  * De Wilde  _et al._ [2006] Y. De Wilde, F. Formanek, R. Carminati, B. Gralak, P.-A. Lemoine, K. Joulain, J.-P. Mulet, Y. Chen, and J.-J. Greffet, Thermal radiation scanning tunnelling microscopy, Nature 444, 740 (2006).
  * Li  _et al._ [2018] C. Li, V. Krachmalnicoff, P. Bouchon, J. Jaeck, N. Bardou, R. Haidar, and Y. De Wilde, Near-field and far-field thermal emission of an individual patch nanoantenna, Physical Review Letters 121, 243901 (2018).
  * Jauffred  _et al._ [2019] L. Jauffred, A. Samadi, H. Klingberg, P. M. Bendix, and L. B. Oddershede, Plasmonic heating of nanostructures, Chemical Reviews 119, 8087 (2019).
  * Cuevas and García-Vidal [2018] J. C. Cuevas and F. J. García-Vidal, Radiative heat transfer, Acs Photonics 5, 3896 (2018).
  * Li and Fan [2018] W. Li and S. Fan, Nanophotonic control of thermal radiation for energy applications, Optics Express 26, 15995 (2018).
  * Makhsiyan  _et al._ [2015] M. Makhsiyan, P. Bouchon, J. Jaeck, J.-L. Pelouard, and R. Haïdar, Shaping the spatial and spectral emissivity at the diffraction limit, Applied Physics Letters 107, 251103 (2015).
  * Baranov  _et al._ [2019] D. G. Baranov, Y. Xiao, I. A. Nechepurenko, A. Krasnok, A. Alù, and M. A. Kats, Nanophotonic engineering of far-field thermal emitters, Nature Materials 18, 920 (2019).
  * Wojszvzyk  _et al._ [2021] L. Wojszvzyk, A. Nguyen, A.-L. Coutrot, C. Zhang, B. Vest, and J.-J. Greffet, An incandescent metasurface for quasimonochromatic polarized mid-wave infrared emission modulated beyond 10 mhz, Nature Communications 12, 1 (2021).
  * García-Esteban  _et al._ [2022] J. J. García-Esteban, J. Bravo-Abad, and J. C. Cuevas, Tunable thermal emission of subwavelength silica ribbons, ACS Photonics 9, 3679 (2022).
  * Nguyen  _et al._ [2022] A. Nguyen, J.-P. Hugonin, E. Garcia-Caurel, A.-L. Coutrot, B. Vest, and J.-J. Greffet, Circularly polarized light emission by incandescent metasurfaces, in  _Plasmonics: Design, Materials, Fabrication, Characterization, and Applications XX_ (SPIE, 2022) p. PC121970I.
  * Cui  _et al._ [2014] Y. Cui, Y. He, Y. Jin, F. Ding, L. Yang, Y. Ye, S. Zhong, Y. Lin, and S. He, Plasmonic and metamaterial structures as electromagnetic absorbers, Laser & Photonics Reviews 8, 495 (2014).
  * Sterl  _et al._ [2021] F. Sterl, E. Herkert, S. Both, T. Weiss, and H. Giessen, Shaping the color and angular appearance of plasmonic metasurfaces with tailored disorder, ACS Nano 15, 10318 (2021).
  * Ogawa and Kimata [2018] S. Ogawa and M. Kimata, Metal-insulator-metal-based plasmonic metamaterial absorbers at visible and infrared wavelengths: a review, Materials 11, 458 (2018).
  * McCool and Fiume [1992] M. McCool and E. Fiume, Hierarchical poisson disk sampling distributions, in  _Graphics interface_ , Vol. 92 (1992) pp. 94–105.
  * Widom [1966] B. Widom, Random sequential addition of hard spheres to a volume, The Journal of Chemical Physics 44, 3888 (1966).
  * Torquato and Stillinger [2003] S. Torquato and F. H. Stillinger, Local density fluctuations, hyperuniformity, and order metrics, Physical Review E 68, 041113 (2003).
  * Vynck  _et al._ [2021] K. Vynck, R. Pierrat, R. Carminati, L. S. Froufe-Pérez, F. Scheffold, R. Sapienza, S. Vignolini, and J. J. Sáenz, Light in correlated disordered media, arXiv preprint arXiv:2106.13892 (2021).
  * Yu  _et al._ [2021] S. Yu, C.-W. Qiu, Y. Chong, S. Torquato, and N. Park, Engineered disorder in photonics, Nature Reviews Materials 6, 226 (2021).
  * Striebel  _et al._ [2017] M. Striebel, J. Wrachtrup, and I. Gerhardt, Absorption and extinction cross sections and photon streamlines in the optical near-field, Scientific Reports 7, 1 (2017).
  * Streed  _et al._ [2012] E. W. Streed, A. Jechow, B. G. Norton, and D. Kielpinski, Absorption imaging of a single atom, Nature Communications 3, 1 (2012).
  * Tretyakov [2014] S. Tretyakov, Maximizing absorption and scattering by dipole particles, Plasmonics 9, 935 (2014).
  * Multiphysics [1998] C. Multiphysics, Introduction to comsol multiphysics®, COMSOL Multiphysics, Burlington, MA, accessed Feb 9, 2018 (1998).


# Data Sketching and Stacking: A Confluence of Two Strategies for Predictive
Inference in Gaussian Process Regressions with High-Dimensional Features

Samuel Gailliot
Department of Statistics, Texas A&M University
and
Rajarshi Guhaniyogi
Department of Statistics, Texas A&M University
and
Roger D. Peng
Department of Statistics and Data Sciences, University of Texas at Austin  The
authors gratefully acknowledge please remember to list all relevant funding
sources in the unblinded version

###### Abstract

This article focuses on drawing computationally-efficient predictive inference
from Gaussian process (GP) regressions with a large number of features when
the response is conditionally independent of the features given the projection
to a noisy low dimensional manifold. Bayesian estimation of the regression
relationship using Markov Chain Monte Carlo and subsequent predictive
inference is computationally prohibitive and may lead to inferential
inaccuracies since accurate variable selection is essentially impossible in
such high-dimensional GP regressions. As an alternative, this article proposes
a strategy to sketch the high-dimensional feature vector with a carefully
constructed sketching matrix, before fitting a GP with the scalar outcome and
the sketched feature vector to draw predictive inference. The analysis is
performed in parallel with many different sketching matrices and smoothing
parameters in different processors, and the predictive inferences are combined
using _Bayesian predictive stacking_. Since posterior predictive distribution
in each processor is analytically tractable, the algorithm allows bypassing
the robustness issues due to convergence and mixing of MCMC chains, leading to
fast implementation with very large number of features. Simulation studies
show superior performance of the proposed approach with a wide variety of
competitors. The approach outperforms competitors in drawing point prediction
with predictive uncertainties of outdoor air pollution from satellite images.

Keywords: Bayesian predictive stacking; feature sketching; Gaussian processes;
high-dimensional features; manifold regression; posterior consistency.

##  1 Introduction

We focus on the problem of drawing predictive inference of a random variable
from a high-dimensional feature vector using “sketching” of the feature vector
when it truly lies on a low-dimensional noisy unknown manifold. In recent
years, there has been a growing literature on “data sketching,” which involves
sketching or compressing the original data before analysis (Halko et al.,
2011, Mahoney et al., 2011, Woodruff et al., 2014, Guhaniyogi and Dunson,
2015, 2016). However, our approach differs from the existing data sketching
literature in two key aspects. Firstly, while most data sketching approaches
aim to reduce the number of data samples, our approach is distinct in that it
maintains the same number of samples but instead reduces the dimensionality of
the feature vector. Secondly, the majority of research in data sketching
focuses on performance evaluation of ordinary and high-dimensional penalized
regression methods with sketched data (Zhang et al., 2013, Dobriban and Liu,
2018, Drineas et al., 2011, Ahfock et al., 2017, Huang, 2018), with only a few
recent articles considering application of data sketching in Bayesian high-
dimensional linear and non-linear regressions (Guhaniyogi and Scheffler, 2021,
Guhaniyogi and Dunson, 2015, 2016). In contrast, our approach leverages the
benefits of data sketching to deliver scalable predictive inference in non-
parametric regressions with a limited sample size and a large number of
features, when the features lie on a noisy low-dimensional manifold.

We consider a regression framework with an outcome
$y\in\mathcal{Y}\subseteq\mathbb{R}$ and a feature vector
${\bm{x}}=(x_{1},...,x_{p})^{T}$ when ${\bm{x}}$ resides on a noisy unknown
manifold, i.e., ${\bm{x}}={\bm{\phi}}({\bm{o}})+{\bm{\eta}}$, where
${\bm{o}}=(o_{1},...,o_{d})^{T}$ is $d$-dimensional co-ordinates for a
manifold $\mathcal{O}\subseteq\mathbb{R}^{p}$,
${\bm{\phi}}(\cdot):\mathbb{R}^{d}\rightarrow\mathbb{R}^{p}$ is a mapping
function such that ${\bm{\phi}}({\bm{o}})\in\mathcal{O}$ and ${\bm{\eta}}$ is
$p$-dimensional noise. Often the complex dependence between $y$ and ${\bm{x}}$
is encoded via co-ordinates of the low-dimensional manifold, i.e.,

$\displaystyle y=h({\bm{o}})+\epsilon,$ (1)

where $h$ is a complex function encoding the true relationship between
response and co-ordinates of the manifold and $\epsilon$ is the error. Since
the manifold $\mathcal{O}$ is unobserved, the co-ordinates ${\bm{o}}$ is
typically unknown. Hence, the common practice is to estimate complex
dependencies between $y$ and ${\bm{x}}$ through a non-linear regression model
given by,

$\displaystyle y=f({\bm{x}})+\epsilon,$ (2)

where $f$ is an unknown regression function and $\epsilon$ is the residual.
When dealing with high-dimensional features, Gaussian process (GP) priors with
an automatic relevance determination (ARD) kernel are commonly used to
estimate the underlying function $f$ with sufficient sparsity assumption in
the relationship between $y$ and ${\bm{x}}$ (Zhao et al., 2018, Jensen et al.,
2021). The estimated $f$ is then employed to predict the response variable.
However, when the number of features reaches the order of a few thousand,
estimation of $f$ with GP-ARD framework is often inaccurate, leading to
unsatisfactory predictive inference. This article proposes an alternative
approach that exclusively focuses on drawing predictive inference on the
response variable $y$, including both point prediction and uncertainty
estimation, using GP regression. We review below a list of existing strategies
to draw predictive inference on $y$ in non-linear regressions before
introducing our approach.

In the literature, a significant line of work follows a two-stage approach for
dealing with high-dimensional features in non-linear manifold regression
tasks. In this approach, the first stage involves constructing a lower-
dimensional representation of the high-dimensional features using manifold
learning techniques. Some popularly employed parsimonious manifold learning
algorithms include Isomap (Tenenbaum et al., 2000), Diffusion Maps (Coifman
and Lafon, 2006), and Laplacian eigenmap (Belkin and Niyogi, 2003). These
algorithms enable the reduction of dimensionality while preserving the
essential characteristics of the data. Additionally, there are model-based
approaches that estimate the unknown Riemannian manifold structure within the
feature vector. These methods utilize techniques such as local PCA (Weingessel
and Hornik, 2000, Arias-Castro et al., 2017) or geometric multiresolution
analysis (Maggioni et al., 2016), and, more recently, spherical basis
functions (Li et al., 2022). Non-linear regression models in the second stage
are based on these projected features in lower-dimensions. However, it is
important to note that such two-stage approaches rely on learning the manifold
structure embedded in the high-dimensional features. While this can be
valuable for understanding the underlying data structure, it adds unnecessary
computational burden when the primary focus is on prediction rather than
inference.

An alternative line of research focuses on estimating the unknown function $f$
using tree-based approaches or deep neural network methods. Tree-based
approaches, such as CART (Denison et al., 1998), BART (Chipman et al., 2010),
and random forest (Breiman, 2001) are based on finding the best splitting
attribute, which can become less efficient as the number of features ($p$)
increases. While there is a growing literature on variable selection within
tree-based methods, such as BART (Bleich et al., 2014) and its variants (Liu
et al., 2021), estimating the true regression function with a large number of
features ($p$ of the order of thousands) can pose challenges. Deep neural
networks with variable selection architecture (Dinh and Ho, 2020) are also not
ideal as they lack predictive uncertainty and struggle to handle the high-
dimensional feature space efficiently.

Bayesian modeling approaches are naturally appealing when the focus is on
quantifying predictive uncertainty. To this end, the more traditional Bayesian
models simultaneously learn the mapping to the lower-dimensional subspace
along with the regression function in the coordinates on this subspace. These
approaches range from Gaussian process latent variable models (GP-LVMs)
(Lawrence and Hyvärinen, 2005, Titsias and Lawrence, 2010) for probabilistic
nonlinear principle component analysis to mixture of factor models (Chen et
al., 2010). However, such methods pose daunting computational challenges with
even moderately large $p$ and sample size due to learning the number and
distribution of latent variables, as well as the mapping functions, while
maintaining identifiability restrictions.

To enhance the time efficiency of the aforementioned approaches, pre-
processing steps are often employed, and two popular pre-processing methods
are feature screening and projection. The feature screening approach
identifies features that exhibit the strongest marginal association with the
response variable. By selecting the features with the highest marginal
association, this approach aims to reduce the dimensionality of the problem.
Feature screening methods are generally straightforward to implement, and it
offers asymptotic guarantees of selecting a superset of important features
(Chen et al., 2018). On the other hand, projection approaches aim to construct
lower-dimensional feature vectors by combining the original high-dimensional
features. One common method in projection approaches is to construct a few
principal components (PCs) from the original $p$-dimensional feature vector.

A naive implementation of the above pre-processing steps is unappealing to our
scenario. For example, in a non-parametric regression with a large number of
correlated features and low signal-to-noise ratio, it may be important to
choose a conservative threshold for screening, which limits the scope of
dimension reduction at this stage. On the other hand, construction of PCs are
agnostic to the relationship between the response and the feature vector.
Instead, we propose an approach that first employs variable screening (Chen et
al., 2018) with a conservative threshold to identify a large subset of
features, typically a few thousand, having the highest non-linear marginal
association with the response. After variable screening, the screened feature
vector is further compressed using a short and fat random sketching matrix
(Mahoney et al., 2011, Drineas et al., 2012). This matrix has a small number
of rows ($m$) and entries that are independently and identically drawn from a
normal distribution. Predictive inference proceeds by fitting a non-parametric
Gaussian process (GP) regression model (Williams and Rasmussen, 2006, Gramacy,
2020) to the scalar outcome and the $m$-dimensional sketched feature vector
after fixing values for the weakly identifiable tuning parameters within the
covariance kernel of the Gaussian processes. The posterior predictive
distribution corresponding to a choice of such tuning parameters and random
sketching matrix comes in a closed form without the need to implement MCMC
sampling, so that one can obtain the predictive distribution extremely rapidly
even in problems with huge numbers of features. To reduce the sensitivity of
predictive inference to the choice of the sketching matrix and the tuning
parameters in GP regression, the model is fit with multiple different choices
of the sketching matrix and parameters. The predictive inferences obtained
from such choices are then aggregated using Bayesian predictive stacking (Yao
et al., 2018) to improve accuracy and robustness.

Stacking is a model aggregation procedure to combine predictions from many
different models (Wolpert, 1992, Breiman, 1996, LeBlanc and Tibshirani, 1996).
In recent years, substantial advancements have been made in Bayesian stacking
methodology, with notable contributions made in Le and Clarke (2017), Yao et
al. (2018), Pavlyshenko (2020), Yao et al. (2022a, b) and the references
therein. However, to the best of our knowledge, the application of stacking in
the context of predictive inference for high-dimensional manifold regression
is currently lacking. While Bayesian model averaging (Raftery et al., 1997) is
most popularly used for aggregating predictive inference from multiple models,
it may be less suited to the stacking procedure in our settings. To see this,
assume that there are $S$ candidate models
$\mathcal{M}=\\{\mathcal{M}_{1},...,\mathcal{M}_{S}\\}$. Bayesian model
comparison typically encounter three different settings: (i)
$\mathcal{M}$-closed where a true data generating model exists and is included
in $\mathcal{M}$; (ii) $\mathcal{M}$-complete where a true model exists but is
not included in $\mathcal{M}$; and (iii) $\mathcal{M}$-open where we do not
assume the existence of a true data generating model. Although Bayesian model
averaging has the advantage of asymptotically identifying the true data
generating model in the first setting, predictive stacking has advantages in
the $\mathcal{M}$-complete and $\mathcal{M}$-open settings. Given that the
true model may not be included in the class of fitted Gaussian process
regression models with randomly sketched features, predictive stacking offers
substantial advantages over model averaging.

In regressions involving high-dimensional features and a large sample size,
Guhaniyogi and Scheffler (2021) propose an approach orthogonal to ours which
exploits random sketching matrices to reduce the sample size rather than the
number of features. A few approaches closely related to ours develop
theoretical bound on predictive accuracy when high-dimensional features are
sketched with random matrices (Guhaniyogi and Dunson, 2015, 2016, Thanei et
al., 2017). These articles tend to include random linear combinations of many
unimportant features, diminishing signal in the analysis, which results in
less than satisfactory predictive performance with massive-dimensional
features. Addressing this issue, Mukhopadhyay and Dunson (2020) proposes novel
constructions of projection matrices tailored to deliver more accurate
predictive inference. These approaches aim to overcome the challenges
associated with sketching high-dimensional features and improve predictive
performance. However, they primarily focus on high-dimensional parametric
regression, and their applicability to non-parametric regression tasks may
require further investigation. Additionally, these approaches address
sensitivity to the choice of sketching matrices by aggregating predictive
inference over many sketching matrices using Bayesian model averaging
technique (Raftery et al., 1997) which is less suitable for prediction than
the stacking approach we employ here, as discussed in the last paragraph.

The rest of the article proceeds as follows. Section 1.1 discusses motivating
dataset on outdoor air pollution and satellite imagery. Section 2 proposes the
model and computational approach for predictive inference in manifold
regression with large number of predictors. Section 3 offers empirical
evaluation of the proposed approach along with its competitors for simulation
studies. Section 4 investigates the proposed approach in drawing predictive
inference of outdoor air pollution concentration from satellite images.
Finally, Section 5 concludes the article with an eye towards future work.

###  1.1 Outdoor Air Pollution Application

As a motivation for the development of our methodology, we consider the
problem of predicting outdoor air pollution concentrations across the United
States. Outdoor air pollution in the U.S. is measured using a network of
ground-based monitors managed by local air quality agencies and the U.S.
Environmental Protection Agency (Environmental Protection Agency, 1996) (EPA).
While the combined network of monitors consists of thousands of locations, the
spatial coverage of the network is actually quite sparse, leaving many areas
of the country without any ground-level data (Apte et al., 2017). Many dense
urban areas only contain one or two monitors, raising a question of whether
such measurements are representative of the burden experienced by all members
of the population. To address the sparsity of the network, there have been
efforts to deploy low-cost sensors across urban areas to fill the gaps. While
such approaches have promise, they are still experimental and ad hoc in
nature, and the sensors themselves can sometimes introduce new measurement
problems (Heffernan et al., 2023).

Remote sensing techniques, which use satellite imagery to predict ground-level
concentrations of outdoor air pollution have the potential to address the
spatial coverage problem because of their constant monitoring of the entire
planet. Traditional approaches have employed aerosol optical depth as a proxy
for such pollutants as fine particulate matter (Paciorek et al., 2008), or
PM2.5. While the previous generation of Earth observation satellites had
excellent spatial coverage, they lacked temporal coverage, typically
revisiting an area of the planet only once every one or two weeks. In
addition, older satellites tended to have lower resolution, making them
difficult to use for predicting air pollution concentrations in dense urban
settings. In recent years, there has been a revolution in the deployment of
satellite constellations, where hundreds of smaller inexpensive satellites
orbit the Earth, providing constant coverage of all areas (Planet Team, 2017).
Furthermore, these satellites have much higher resolution, allowing for more
detailed examination of areas of interest.

Given the recent emergence of data from satellite constellations, there is
still a question of how best to use them for the purpose of predicting ground-
level air pollution concentrations. For this application, we focus on
predicting fine particulate matter pollution (PM2.5) from multi-band satellite
images. For ground-truth information we use the EPA’s network of monitors to
provide valid PM2.5 measurements. The combination of high-resolution spatial
and temporal coverage of the entire U.S. with novel statistical prediction
approaches has the potential to dramatically increase the monitoring of
outdoor air pollution and its subsequent health effects.

##  2 Our Approach: Stacked Gaussian Process Regression

Let $\mathcal{D}_{n}=\\{({\bm{x}}_{i}^{T},y_{i}):i=1,...,n\\}$ be a dataset
containing $n$ observations each with a $p$-variate feature
${\bm{x}}_{i}=(x_{i,1},...,x_{i,p})^{T}$ and a scalar-valued response $y_{i}$.
We assume $n$ is moderately large and $p$ is large. The feature vector
${\bm{x}}_{i}$ lies on an unknown noisy manifold
$\mathcal{O}\subseteq\mathbb{R}^{p}$ with $d$-dimensional latent co-ordinates
${\bm{o}}_{i}$ (i.e.,
${\bm{x}}_{i}={\bm{\phi}}({\bm{o}}_{i})+{\bm{\eta}}_{i}$,
${\bm{\phi}}({\bm{o}}_{i})\in\mathcal{O}$). We assume a nonlinear regression
relationship between $y_{i}$ and ${\bm{x}}_{i}$, and approximate the density
of $y_{i}$ by sketching the high-dimensional feature vector ${\bm{x}}_{i}$ to
lower-dimensions using a sketching matrix ${\bm{P}}_{n}$ as follows

$\displaystyle
y_{i}=f({\bm{P}}_{n}{\bm{x}}_{i})+\epsilon_{i},\>\>\epsilon_{i}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}N(0,\xi^{2}),$ (3)

with $\xi^{2}$ as the noise variance and $f(\cdot)$ as an unknown continuous
function in the Holder class of smoothness $s$. Discussion on the choice of
the sketching matrix ${\bm{P}}_{n}\in\mathbb{R}^{m\times p}$ is provided in
Section 2.1.

###  2.1 Choice of the Sketching Matrix

The sketching matrix ${\bm{P}}_{n}\in\mathbb{R}^{m\times p}$ embeds
$p$-dimensional features ${\bm{x}}_{i}$ into $m$ dimensions while not throwing
away excessive amounts of information. The most popular linear embedding is
obtained from the singular value decomposition (SVD) of
${\bm{X}}=[{\bm{x}}_{1}:\cdots:{\bm{x}}_{n}]^{T}$, but they are problematic to
estimate when $p>>n$. In contrast, random sketching matrices are often used to
embed the high-dimensional features to a random subspace, and appropriate
choices of the random matrices allow distances between samples to be
approximately preserved (Li and Gu, 2017).

The direct application of sketching matrices on high-dimensional features is
unappealing, as it constructs random linear combinations of many unimportant
features, diminishing the signal in the analysis. As an alternative approach,
we design a sketching matrix that constructs random linear combinations of
features with the highest marginal association with the response. To identify
these features, we perform nonparametric B-spline regression of $y_{i}$ onto
each component of ${\bm{x}}_{i}$ separately. The order of importance of the
features is determined, in descending order, by the residual sum of squares of
the marginal nonparametric regressions. Features with a residual sum of
squares greater than a user-defined threshold are considered important
features related to the response. We adopt a conservative threshold following
Fan et al. (2014) to select a large superset of important features, which
allows for joint contributions of features in explaining the response.

Let $\mathcal{I}$ correspond to the indices of the features chosen with
marginal screening and $\bar{\mathcal{I}}$ be the indices of the features
screened out through this procedure, such that
$\mathcal{I}\cup\bar{\mathcal{I}}=\\{1,...,p\\}$. Let ${\bm{E}}_{n}$ be a
permutation matrix such that
${\bm{E}}_{n}{\bm{x}}=({\bm{x}}_{\mathcal{I}}^{T},{\bm{x}}_{\bar{\mathcal{I}}}^{T})^{T}$.
We construct a matrix ${\bm{R}}_{n}=[{\bm{R}}_{n,1}:{\bm{R}}_{n,2}]$ where
${\bm{R}}_{n,2}={\bm{0}}_{m\times(p-|\mathcal{I}|)}$ and ${\bm{R}}_{n,1}$ is
an $m\times|\mathcal{I}|$ matrix with entries drawn independently from N(0,1),
following the literature on random sketching matrices (Baraniuk et al., 2008).
The resulting sketching matrix is given by
${\bm{P}}_{n}={\bm{R}}_{n}{\bm{E}}_{n}$.

###  2.2 Prior, Posterior and Posterior Predictive Distributions

Following a Bayesian approach, we assign a zero-centered Gaussian process
prior on the unknown regression function $f(\cdot)$, denoted by $f(\cdot)\sim
GP(0,\sigma^{2}\delta_{\theta})$. Here, $\delta_{\theta}$ corresponds to an
exponential correlation kernel
$\delta_{\theta}({\bm{x}}_{i},{\bm{x}}_{j})=\exp(-\theta||{\bm{x}}_{i}-{\bm{x}}_{j}||)$
involving the length-scale parameter $\theta$. The parameter $\sigma^{2}$ is
the signal variance parameter and $||\cdot||$ denotes the Euclidean norm. A
significant finding by Yang and Dunson (2016) establishes that when features
${\bm{x}}_{i}$ lie on a $d$-dimensional manifold $\mathcal{O}$, the minimax
optimal rate of $n^{-2s/(2s+d)}$ (adapted to the dimension of the manifold)
can be achieved in estimating $f$ through an appropriate choice of prior
distributions on $\theta$ and $\sigma^{2}$. However, in practical scenarios,
features may not exactly lie on a manifold due to noise and data corruption,
as assumed in our setting. In such instances, the application of random
compression, denoted as ${\bm{P}}_{n}{\bm{x}}_{i}$, aids in denoising the
features. The de-noised compressed features ${\bm{P}}_{n}{\bm{x}}_{i}$ exhibit
a higher concentration around the manifold compared to the original features
${\bm{x}}_{i}$. With this enhanced concentration, the theory presented by Yang
and Dunson (2016) suggests that an appropriate GP prior can yield excellent
performance. In addition to denoising, the compression of the high-dimensional
feature vector has a major advantage in avoiding the estimation of a geodesic
distance along the unknown manifold $\mathcal{O}$ between any two feature
vectors ${\bm{x}}_{i}$ and ${\bm{x}}_{i^{\prime}}$.

In practical applications, utilizing the recommended prior distributions on
$\theta$ and $\sigma^{2}$ from Yang and Dunson (2016) requires computationally
expensive Markov Chain Monte Carlo (MCMC) sampling, mainly due to the weak
identifiability of $\theta$. The posterior computation of $\theta$ typically
entails meticulous tuning, especially when dealing with high-dimensional
features, imposing a significant computational burden. This article introduces
an alternative approach that enables exact predictive inference from the
model, entirely bypassing the MCMC algorithm in model estimation. The details
of the strategy are elaborated below.

Denote
${\bm{f}}=(f({\bm{P}}_{n}{\bm{x}}_{1}),...,f({\bm{P}}_{n}{\bm{x}}_{n}))^{T}$
as the vector consisting of the function $f$ evaluated at the sketched
features ${\bm{P}}_{n}{\bm{x}}_{1}$,…,${\bm{P}}_{n}{\bm{x}}_{n}$ and
${\bm{C}}$ as an $n\times n$ covariance matrix with $(i,j)$th entry
$\delta_{\theta}({\bm{P}}_{n}{\bm{x}}_{i},{\bm{P}}_{n}{\bm{x}}_{j})$. With
${\bm{y}}=(y_{1},...,y_{n})^{T}$ as the response vector, a customary Bayesian
hierarchical model is constructed as

$\displaystyle{\bm{y}}|{\bm{f}},\xi^{2}\sim
N({\bm{f}},\xi^{2}{\bm{I}}),\>\>\>({\bm{f}}|\xi^{2})\sim
N(0,\xi^{2}\psi^{2}{\bm{C}}),\>\>\>\>\pi(\xi^{2})\propto\frac{1}{\xi^{2}},$

where we fix the length-scale parameter $\theta$ and the signal-to-noise
variance ratio $\psi^{2}=\frac{\sigma^{2}}{\xi^{2}}$. This ensures closed-form
conjugate marginal posterior and posterior predictive distributions. More
specifically, the marginal posterior distribution of $\xi^{2}$, given the
projection matrix ${\bm{P}}_{n}$, $\theta$, $\psi^{2}$ and $\mathcal{D}_{n}$,
is inverse gamma with parameters $a=n/2$ and
$b={\bm{y}}^{T}(\psi^{2}{\bm{C}}+{\bm{I}})^{-1}{\bm{y}}/2$. The marginal
posterior distribution of ${\bm{f}}$, given ${\bm{P}}_{n}$, $\psi^{2}$,
$\theta$ and $\mathcal{D}_{n}$, follows a scaled $n$-variate $t$ distribution
with degrees of freedom $n$, location ${\bm{\mu}}_{t}$ and scale matrix
${\bm{\Sigma}}_{t}$, denoted by $t_{n}({\bm{\mu}}_{t},{\bm{\Sigma}}_{t})$,
where
${\bm{\mu}}_{t}=({\bm{I}}+{\bm{C}}^{-1}/\psi^{2})^{-1}{\bm{y}},\>\>\>{\bm{\Sigma}}_{t}=(2b/n)({\bm{I}}+{\bm{C}}^{-1}/\psi^{2})^{-1}.$
Consider prediction for the response at $n_{new}$ data points with
corresponding covariates
$\tilde{{\bm{x}}_{1}},...,\tilde{{\bm{x}}}_{n_{new}}$. Let
${\bm{C}}_{new,new}$ and ${\bm{C}}_{new,old}$ denote $n_{new}\times n_{new}$
and $n_{new}\times n$ matrices with $(i,j)$th elements
$\delta_{\theta}({\bm{P}}_{n}\tilde{{\bm{x}}}_{i},{\bm{P}}_{n}\tilde{{\bm{x}}}_{j})$
and
$\delta_{\theta}({\bm{P}}_{n}\tilde{{\bm{x}}}_{i},{\bm{P}}_{n}{\bm{x}}_{j})$,
respectively. The posterior predictive distribution of the response
$\tilde{{\bm{y}}}_{new}=(\tilde{y}_{1},...,\tilde{y}_{n_{new}})^{T}$ given
$\tilde{{\bm{x}}_{1}},...,\tilde{{\bm{x}}}_{n_{new}}$, ${\bm{P}}_{n}$,
$\theta$, $\psi^{2}$ and $\mathcal{D}_{n}$, marginalizing out
$({\bm{f}},\xi^{2})$, follows a scaled $n_{new}$-variate t-distribution
$t_{n_{new}}(\tilde{{\bm{\mu}}}_{t},\tilde{{\bm{\Sigma}}}_{t})$, where

$\displaystyle\tilde{{\bm{\mu}}}_{t}$
$\displaystyle=\psi^{2}{\bm{C}}_{new,old}({\bm{I}}+\psi^{2}{\bm{C}})^{-1}{\bm{y}}$
$\displaystyle\tilde{{\bm{\Sigma}}}_{t}$
$\displaystyle=(2b/n)\left[{\bm{I}}+\psi^{2}{\bm{C}}_{new,new}-\psi^{4}{\bm{C}}_{new,old}({\bm{I}}+\psi^{2}{\bm{C}})^{-1}{\bm{C}}_{new,old}^{T}\right].$
(4)

Since the posterior predictive distribution is available in closed form,
Bayesian inference can proceed from exact posterior samples.

This tractability is only possible if the length-scale parameter $\theta$ and
the signal-to-noise variance ratio $\psi^{2}$ are fixed. While it is possible
to estimate their full posterior distributions through expensive Markov Chain
Monte Carlo (MCMC) sampling, these parameters are inconsistently estimable for
the general Matern class of correlation functions (Zhang, 2004) often
resulting in poorer convergence. Therefore, for the chosen sketching matrix
${\bm{P}}_{n}$, we obtain $(\theta,\psi^{2})$ such that

$\displaystyle\mbox{max}_{\theta,\psi^{2}}f(\theta,\psi^{2}|{\bm{P}}_{n},{\bm{y}})$
$\displaystyle\propto\mbox{max}_{\theta,\psi^{2}}\frac{1}{|\psi^{2}{\bm{C}}+{\bm{I}}|^{\frac{1}{2}}}\frac{2^{\frac{n}{2}}\Gamma(\frac{n}{2})}{\left[{\bm{y}}^{\prime}(\psi^{2}{\bm{C}}+{\bm{I}})^{-1}{\bm{y}}\right]^{\frac{n}{2}}(\sqrt{2\pi})^{n}}$
(5)

Our approach will conduct exact predictive inference using the closed form
predictive distribution in (2.2) and stack the predictive inference over
different fixed values of $\\{{\bm{P}}_{n},\theta,\psi^{2}\\}$.

####  2.2.1 Stacking of Predictive Distributions

Let $\mathcal{M}_{k}$ represent the fitted model (3) with
${{\bm{P}}_{n}^{(k)},\theta^{(k)},\psi^{2(k)}}$ for $k=1,...,K$. While the
sketching matrix ${\bm{P}}_{n}^{(k)}$ is randomly generated for each $k$,
$\theta^{(k)}$ and $\psi^{2(k)}$ are obtained following equation (5) for the
choice of ${\bm{P}}_{n}^{(k)}$. Employing the generalized Bayesian stacking
framework proposed by Yao et al. (2018), we implement a stacking procedure
over the predictive distribution obtained from each $\mathcal{M}_{k}$. Let
$p(\tilde{{\bm{y}}}_{new}|{\bm{y}},\mathcal{M}_{k})$ denote the predictive
distribution under model $\mathcal{M}_{k}$, and
$p_{t}(\tilde{{\bm{y}}}_{new}|{\bm{y}})$ denote the true predictive
distribution. Our objective is to determine the distribution in the convex
hull
$\mathcal{C}=\\{\sum_{k=1}^{K}w_{k}p(\cdot|{\bm{y}},\mathcal{M}_{k}):w_{k}\in\mathcal{S}_{1}^{K}\\}$,
where $\mathcal{S}_{1}^{K}=\\{{\bm{w}}\in[0,1]^{K}:\sum_{k=1}^{K}w_{k}=1\\}$,
that is optimal with respect to some proper scoring function. Using the
logarithmic score, which corresponds to the KL divergence, we seek to find the
vector $\tilde{{\bm{w}}}=(\tilde{w}_{1},\ldots,\tilde{w}_{K})$ such that

$\tilde{{\bm{w}}}=\max_{{\bm{w}}\in\mathcal{S}_{1}^{K}}\frac{1}{n}\sum_{i=1}^{n}\log\left(\sum_{k=1}^{K}w_{k}p_{k,-i}(y_{i})\right),$
(6)

where ${\bm{y}}_{-i}=(y_{j}:j\neq i,\>j=1,...,n)^{T}$ and
$p_{k,-i}(y_{i})=p(y_{i}|{\bm{y}}_{-i},\mathcal{M}_{k})$ has a closed from
univariate $t$-distribution with parameters $\tilde{\mu}_{-i}^{(k)}$ and
$\tilde{\Sigma}_{-i}^{(k)}$ obtained using the formula for posterior
predictive distribution given in equation (2.2). In practice, calculating the
predictive densities $p_{k,-i}(y_{i})$ one at a time is computationally
expensive as the calculation of $\tilde{\Sigma}_{-i}^{(k)}$ requires inverting
an $(n-1)\times(n-1)$ matrix for every $k=1,...,K$ and $i=1,...,n$. To avoid
this, we randomly split the data into $S=10$ disjoint folds of approximately
equal size,
$({\bm{y}}_{(1)},{\bm{X}}_{(1)}),\ldots,({\bm{y}}_{(S)},{\bm{X}}_{(S)})$, and
compute
${\bm{y}}_{(s)}|{\bm{y}}_{(1)},...,{\bm{y}}_{(s-1)},{\bm{y}}_{(s+1)},...,{\bm{y}}_{(S)}$
for every $s=1,...,S$, which follows a multivariate t-distribution with
parameters obtained using equation (2.2). If the $i$th sample belongs to the
$s$th fold, we will replace $p_{k,-i}(y_{i})$ in (6) by $p_{k,(s)}(y_{i})$,
where $p_{k,(s)}(y_{i})$ represents the marginal distribution of $y_{i}$ from
${\bm{y}}_{(s)}|{\bm{y}}_{(1)},...,{\bm{y}}_{(s-1)},{\bm{y}}_{(s+1)},...,{\bm{y}}_{(S)}$.
This strategy requires inverting an $(n-n/S)\times(n-n/S)$ matrix only $S$
times (assuming that all folds are of equal size), which leads to substantial
computational benefits. No analytical solution to this non-convex constrained
optimization problem in (6) is available, but first and second derivatives are
easily obtained to construct an iterative optimizer. The optimal distribution
provides a pseudo posterior predictive distribution given by
$\tilde{p}(\tilde{{\bm{y}}}_{new}|{\bm{y}})=\sum_{k=1}^{K}\tilde{w}_{k}t_{n_{new}}(\tilde{{\bm{\mu}}}_{t}^{(k)},\tilde{{\bm{\Sigma}}}_{t}^{(k)})$,
where $\tilde{{\bm{\mu}}}_{t}^{(k)}$ and $\tilde{{\bm{\Sigma}}}_{t}^{(k)}$ are
obtained from equation (2.2) by evaluating $\tilde{{\bm{\mu}}}_{t}$ and
$\tilde{{\bm{\Sigma}}}_{t}$ at ${\bm{P}}_{n}^{(k)},\theta^{(k)},\psi^{(k)}$.
The pseudo posterior predictive distribution is further used to draw point
prediction and 95% predictive interval to quantify predictive uncertainty.
Figure 1 offers a flowchart outlining the proposed framework.

Figure 1: Flowchart representing the Sketched Gaussian process regression
framework.

While Bayesian model averaging (BMA) is a common method for combining multiple
distributions, its applicability in our context is limited for several
reasons. Firstly, the fitted models (3) with randomly sketched features are
likely to deviate from the true model, placing us outside the
$\mathcal{M}$-closed setting where BMA is optimal. Additionally, stacking is
designed to determine weights for optimal prediction, whereas asymptotically,
BMA assigns full weight to the “best” single model closest in KL divergence to
the true model (Yao et al., 2022a). However, when the true model lies outside
the space of fitted models, it may be more advantageous to leverage multiple
models in predictive inference. Subsequent empirical experiments demonstrate
stacking as a powerful tool for drawing posterior predictive inference in our
setting.

##  3 Simulation Study

We evaluate the performance of the proposed Sketched Gaussian Process (SkGP)
regression across various simulation scenarios, exploring different structure
of the manifold ($\mathcal{O}$), different feature dimensions ($p$) and noise
levels in the features ($\tau^{2}$) to analyze their impact. In all
simulations, the out-of-sample predictive performance of the proposed SkGP
regression is compared with that of uncompressed Gaussian Process (GP),
Bayesian Additive Regression Trees (BART) (Chipman et al., 2010), Random
Forests (RF) (Breiman, 2001), and deep neural network (NN). We also explore
sketched versions of BART and RF, referred to as Sketched BART (SkBART) and
Sketched Random Forest (SkRF), respectively, where a single projection matrix
is generated to sketch the features, allowing for faster implementation. Each
of these methods are applied on $|\mathcal{I}|=1000$ screened features having
highest marginal association with the response. As a default in this analysis,
we set $m=$ 60\. We offer detailed sensitivity analysis with varying choices
of the number of screened features $|\mathcal{I}|$ and the sketching
dimensions $m$.

###  3.1 Simulated Data Generation

During data simulation, we explore specific scenarios where the response
distribution follows a nonlinear function of $d$-dimensional coordinates for a
manifold $\mathcal{O}\subseteq\mathbb{R}^{p}$, embedded in a high-dimensional
ambient space. Two distinct choices for $\mathcal{O}$ and their corresponding
response distributions are simulated.
$\mathcal{O}$ is a swiss roll and $d=2$. For the swiss roll, we sample
manifold coordinates, $o_{1}\sim U(\frac{3\pi}{2},\frac{9\pi}{2})$, $o_{2}\sim
U(0,3)$. A high dimensional feature ${\bm{x}}=(x_{1},\ldots,x_{p})$ is
simulated according to
$x_{1}=o_{1}\cos(o_{1})+\eta_{1},\;x_{2}=o_{2}+\eta_{2},\;x_{3}=o_{1}\sin(o_{1})+\eta_{3},\;x_{i}=\eta_{i},\;i\geq
4$. The response $y$ have a non-linear relationship with these features and is
simulated following,

$\displaystyle y=sin(5\pi o_{1})+o_{2}^{2}+\epsilon,\>\>\epsilon\sim
N(0,0.02^{2}),$ (7)

where $\eta_{1},\ldots,\eta_{p}\sim N(0,\tau^{2})$. Notably, ${\bm{x}}$ and
$y$ are conditionally independent given $o_{1},o_{2}$ which is the low-
dimensional signal manifold. In particular, ${\bm{x}}$ lives on a (noise
corrupted) swiss roll embedded in a $p$-dimensional ambient space (see Figure
2(a)), but y is only a function of coordinates along the swiss roll
$\mathcal{O}$.
$\mathcal{O}$ is a torus and $d=3$. For the torus, we consider
$x_{1}=o_{1}+\eta_{1},x_{2}=o_{2}+\eta_{2}$ and $x_{3}=o_{3}+\eta_{3}$ where
$o_{1},o_{2},o_{3}$ lie on a three dimensional torus with interior radius $1$
and exterior radius $3$ (see Figure 2(b)), such that
$(3-\sqrt{o_{1}^{2}+o_{2}^{2}})^{2}+o_{3}^{2}=1$, and set $x_{i}=\eta_{i}$ for
$i\geq 4$. The feature noise $\eta_{1},...,\eta_{n}$ are generated i.i.d. from
$N(0,\tau^{2})$. The response is generated as,

$\displaystyle y=o_{2}^{2}+\sin(5\pi o_{3})+\epsilon,\;\epsilon\sim
N(0,0.1^{2}).$

The geodesic distance between two points on both a swiss roll and a torus can
substantially differ from their Euclidean distance in the ambient space
$\mathbb{R}^{p}$. The swiss roll, in particular, poses a challenging setup for
SkGP, as points on $\mathcal{O}$ that are close in a Euclidean sense can be
quite far in a geodesic sense.

To assess the impact of the number of features ($p$) and noise levels of the
features ($\tau^{2}$) on the performance of the competitors, various
simulation scenarios are considered by varying $p=2000,10000$ and
$\tau^{2}=0.01,0.03,0.05,0.1$. For each of these simulation scenarios, 50
datasets are generated, and metrics such as mean squared prediction error
(MSPE), coverage, and lengths of 95% predictive intervals (PI) are calculated
across all replicates. All simulations set the sample size $n=100$ and the
number of predicted samples $n_{new}=100$.

(a) (b)

Figure 2: Manifolds embedded in noisy ambient dimensions

###  3.2 Point prediction

Tables 1 and 2 display the MSPE averaged over $50$ replications for all the
competing methods in the swiss roll and torus examples, respectively. Values
in parentheses represent the standard error of MSPE over $50$ replicates.

Both Tables 1 and 2 show that incorporating randomly sketched features into
the GP model within the SkGP framework yields strong predictive performance,
significantly surpassing the performance of the neural network. For both
$p=2000$ and $p=10000$, when the manifold is affected by low noise, SkGP
significantly outperforms GP, BART, and RF with unsketched features. While
SkBART emerges as the second-best performer in scenarios with very low noise
in the manifold ($\tau^{2}=0.01$), its performance declines notably with an
increase in the noise level in the features. In comparison, SkGP effectively
mitigates the impact of noise in the features, but there exists a tipping
point (depending on the structure of the underlying manifold $\mathcal{O}$ and
sample size $n$) where noise distorts the manifold excessively, causing SkGP
to perform similarly to other competitors. This is observed in the MSPE values
corresponding to $\tau^{2}=0.1$. Among the sketched competitors, SkRF exhibits
notably inferior performance compared to both SkGP and SkBART in all
simulation examples. While theoretically the performance of SkGP should remain
similar for both $p=2000$ and $p=10000$ when all features lie exactly on a
low-dimensional manifold, in practice, we observe a significant decline in the
performance of SkGP with an increase in $p$. This is attributed to the
substantial impact of noise corruption on the manifold, influencing predictive
performance.

Swiss Roll |  | Noise |  |  |
---|---|---|---|---|---
|  | 0.01 | 0.03 | 0.05 | 0.1
| SkGP | 0.96 (0.39) | 1.27 (0.53) | 1.74 (0.431) | 3.26 (0.61)
| GP | 1.28 (1.03) | 1.75 (1.31) | 2.49 (1.12) | 4.81 (0.85)
| BART | 2.31 (0.57) | 2.28 (0.49) | 2.35 (0.49) | 2.57 (0.66)
p = 2000 | SkBART | 0.91 (0.37) | 1.62 (0.51) | 2.63 (0.90) | 5.17 (1.07)
| RF | 6.92 (0.84) | 6.87 (0.87) | 6.92 (0.91) | 6.97 (0.86)
| SkRF | 0.99 (0.56) | 2.05 (0.85) | 3.28 (0.91) | 5.84 (0.97)
| NN | 3.52 (0.61) | 6.66 (0.93) | 7.41 (1.13) | 8.47 (1.08)
| SkGP | 1.64 (0.48) | 2.55 (0.48) | 3.57 (0.57) | 4.54 (0.91)
| GP | 5.19 (1.00) | 5.65 (1.00) | 6.16 (0.99) | 7.08 (0.84)
| BART | 4.17 (1.12) | 4.07 (0.85) | 4.33 (1.06) | 4.37 (0.84)
p = 10,000 | SkBART | 2.38 (0.99) | 5.33 (0.96) | 6.34 (0.86) | 7.31 (0.87)
| RF | 7.32 (0.84) | 7.30 (0.89) | 7.32 (0.95) | 7.35 (0.88)
| SkRF | 3.57 (0.86) | 5.72 (0.96) | 6.66 (0.90) | 7.46 (0.88)
| NN | 7.32 (1.17) | 8.88 (1.58) | 10.15 (1.43) | 10.80 (1.77)
Table 1: Averaged Mean squared Prediction Error (MSPE) over $50$ replications are shown for the competing models in swiss roll example. Standard errors are presented within parenthesis. Torus |  | Noise |  |  |
---|---|---|---|---|---
|  | 0.01 | 0.03 | 0.05 | 0.1
| SkGP | 0.153 (0.036) | 0.281 (0.072) | 0.426 (0.108) | 0.884 (0.127)
| GP | 0.210 (0.081) | 0.308 (0.085) | 0.468 (0.104) | 0.817 (0.119)
| BART | 0.932 (0.132) | 0.992 (0.137) | 0.970 (0.137) | 0.916 (0.152)
p = 2000 | SkBART | 0.201 (0.084) | 0.368 (0.092) | 0.645 (0.151) | 0.990 (0.139)
| RF | 0.957 (0.132) | 1.01 (0.124) | 0.984 (0.124) | 0.937 (0.132)
| SkRF | 0.210 (0.091) | 0.417 (0.103) | 0.641 (0.132) | 0.896 (0.127)
| NN | 0.547 (0.139) | 0.881 (0.122) | 0.977 (0.086) | 1.082 (0.130)
| SkGP | 0.196 (0.048) | 0.503 (0.103) | 0.908 (0.121) | 0.968 (0.136)
| GP | 0.237 (0.077) | 0.649 (0.110) | 0.991 (0.129) | 0.956 (0.132)
| BART | 0.981 (0.129) | 1.04 (0.128) | 1.01 (0.129) | 0.958 (0.132)
p = 10,000 | SkBART | 0.228 (0.096) | 0.897 (0.152) | 1.01 (0.130) | 0.995 (0.136)
| RF | 0.976 (0.136) | 1.03 (0.128) | 1.00 (0.125) | 0.953 (0.129)
| SkRF | 0.328 (0.102) | 0.79 (0.139) | 0.941 (0.131) | 0.975 (0.131)
| NN | 0.919 (0.114) | 1.010 (0.117) | 1.053 (0.125) | 1.094 (0.118)
Table 2: Averaged Mean squared Prediction Error (MSPE) over $50$ replications
are shown for the competing models in torus example. Standard errors are
presented within parenthesis.

###  3.3 Predictive Uncertainty

To evaluate quality of predictive uncertainty, we calculate the coverage and
length of 95% predictive intervals (PI) for SkGP and other competitors. While
frequentist methods, like SkRF and RF, do not inherently provide coverage
probabilities with point estimates, we employ a two-stage plug-in approach for
them: (i) estimate the regression function in the first stage, and (ii)
construct 95% PI based on the normal distribution centered on the predictive
mean from the regression model, with variance equal to the estimated variance
in the residuals. Coverage probability boxplots over $50$ replications for all
simulation cases in the swiss roll example and torus example are presented in
Figures 3 and 4, respectively. Figure 5 displays the median lengths of the 95%
PI for all competitors for all simulation cases in both the swiss roll and
torus examples.

The results indicate that in all simulation scenarios, the coverage of 95%
predictive intervals (PI) for SkGP is close to the nominal level. Although the
intervals tend to widen with both increasing noise in the manifold (i.e.,
higher $\tau^{2}$) and an increase in the number of features $p$, the effect
is less pronounced with $p$ compared to $\tau^{2}$. Both BART and SkBART
exhibit poor coverage, with significantly narrower PIs. RF and SkRF show
undercoverage (around 80% coverage) and wider PIs than SkGP. In the swiss roll
example, the coverage of 95% PI for GP is similar to that of SkGP, but the
intervals from GP are approximately twice as wide as those from SkGP. In the
torus example, GP and SkGP perform similarly for $p=2000$. However, when
$p=10000$, SkGP demonstrates much narrower predictive intervals than GP, with
a similar coverage, as the noise increases in the manifold. Overall, the
results suggest that SkGP is more precise in terms of predictive uncertainty
and robust compared to its competitors concerning the noise in the manifold.

Swiss Roll: $p=2000$

(a) $\tau^{2}=0.01$

(b) $\tau^{2}=0.03$

(c) $\tau^{2}=0.05$

(d) $\tau^{2}=0.1$

Swiss Roll: $p=10,000$

(e) $\tau^{2}=0.01$

(f) $\tau^{2}=0.03$

(g) $\tau^{2}=0.05$

(h) $\tau^{2}=0.1$

Figure 3: Coverage of 95% predictive interval for Swiss Roll Simulations

Torus: $p=2000$

(a) $\tau^{2}=0.01$

(b) $\tau^{2}=0.03$

(c) $\tau^{2}=0.05$

(d) $\tau^{2}=0.1$

Torus: $p=10,000$

(e) $\tau^{2}=0.01$

(f) $\tau^{2}=0.03$

(g) $\tau^{2}=0.05$

(h) $\tau^{2}=0.1$

Figure 4: Coverage of 95% predictive interval for Torus Simulations

(a) Swiss roll (b) Torus

Figure 5: Length of 95% predictive intervals for all competitors in all
simulation settings.

###  3.4 Computation Time

The main objective in developing SkGP was to improve computational scalability
in large $p$ settings. For a specific choice of ${\bm{P}}_{n},\theta$ and
$\psi^{2}$, the computation time for SkGP is primarily influenced by two
factors: (a) computing the inverse of an $n\times n$ matrix; and (b)
multiplying an $m\times|\mathcal{I}|$ matrix with an $|\mathcal{I}|\times n$
matrix. Steps (a) and (b) entail computational complexities of order $n^{3}$
and $mn|\mathcal{I}|$, respectively. Since the posterior predictive
distribution is available in closed forms, these computations are only needed
once for a specific choice of ${\bm{P}}_{n},\theta$ and $\psi^{2}$. We will
parallelize the computation across various choices of
${\bm{P}}_{n},\theta,\psi^{2}$ on different CPUs. The combination step using
stacking requires inverting $S$ matrices each of dimension
$(n-n/S)\times(n-n/S)$, incurring a complexity of the order $S(n-n/S)^{3}$.
Since the focus of this article is on moderate $n$, all computational steps
are extremely efficient leading to rapid computation of SkGP.

Figure 6(a) shows the computation times when the number of features increase
and sketching dimension is held fixed ($m=60$). Computation times for non-
sketched tree based methods increase linearly with the number of features,
while computation time for sketched methods remain constant, as they only
depend on the number of screened features $\mathcal{I}$. Figure 6(b) shows the
computation times as the sketching dimension is increased while the original
number of screened features is held constant ($p=1000$). Considering that the
non-sketched tree based methods are not dependent on sketching dimension $m$,
their computation times remain constant. The computation times for the
sketched methods increase linearly with sketching dimension. Importantly, SkGP
achieves computation time comparable to frequentist approaches, yet being able
to allow principled Bayesian predictive inference.

(a) Vary $p$, fix $m=60$ (b) Vary $m$, fix $|\mathcal{I}|=1000$

Figure 6: The left panel shows the computation time for competitors by fixing
the sketching dimension ($m$), while varying the number of features. The right
panel shows the computation time for competitors by fixing the number of
screened features $|\mathcal{I}|=1000$, while varying $m.$

###  3.5 Sensitivity to the choice of $m$ and $|\mathcal{I}|$

We present investigation into the choice of the number of features
$|\mathcal{I}|$ included through highest marginal association with the
response and the dimension of the sketching matrix $m$ applied to this
$|\mathcal{I}|$-dimensional feature vector. Figure 7 illustrates the MSPE
values for the swiss roll example with varying numbers of included features.
Considering the small true dimensions of the swiss roll manifold and the fact
that the response is related to ${\bm{x}}$ only through the manifold in the
swiss roll example, the inclusion of more redundant features in the regression
leads to a performance loss, as evidenced by the increasing MSPE values.
However, this decline in performance is more pronounced when the swiss roll is
affected by noise with higher variance. This aligns with the fact that the
accuracy of estimating the regression function depends solely on the intrinsic
dimension of the swiss roll and is unaffected by the number of screened
features when the features lie on a manifold. When the noise variance is low,
resulting in features that approximately lie on the swiss roll, the
performance does not change significantly with variations in the number of
screened features.

(a)

(b)

(c)

(d)

Figure 7: We show the number of included features through the marginal
association analysis vis-a-vis predictive accuracy. The plots are presented
for the swiss roll example. Figure 8: Distribution of MSPE as sketching
dimension increases

Figure 8 illustrates the impact of the sketching dimension $m$ on the
performance of SkGP. The figure indicates that as the sketching dimension $m$
increases, MSPE decreases up to a certain point, contingent on the sample size
and the structure of the manifold, after which it starts increasing again.
This observation is reasonable as, with a low-dimensional manifold
$\mathcal{O}$, increasing the sketching dimension $m$ introduces redundant
randomly sketched features to the model, leading to a natural decline in
performance. In practice, we observed that a sketching dimension of around
$m\sim 50$ works well for a diverse range of simulation examples when the
intrinsic dimensionality $d$ of the manifold is low.

##  4 Analysis of Outdoor Air Pollution Data with Satellite Images

We will apply the proposed approach and compare it with relevant competitors
in the analysis of air quality using multi-band satellite images over time
(see Section 1.1). The air quality dataset consists of measurements taken at
the EPA federal reference monitor in Las Vegas, Nevada, spanning almost daily
measurements, sometimes with multiple readings in a day, from January 2019 to
July 2022. This results in $1667$ air quality measurements, as depicted in the
first row of Figure 9. For each air quality sample, multi-band satellite
images covering the location of the air quality monitor have been acquired,
including four wavelength bands: blue, green, red, and near-infrared. The left
panel of Figure 10 displays a near-infrared image on a representative day.
These data were obtained from Planet using version 1 of their PlanetScope
instrument (Planet Team, 2017). Notably, multi-band satellite imagery data are
easily obtainable, whereas the installation of monitors measuring air quality
is expensive. Hence, a key scientific goal is to predict air quality readings
given the high-dimensional multi-band images. To achieve this, at each time
point, the $128\times 128=16384$ pixels of the four bands of the multi-band
images are vectorized and concatenated into a $65536$ dimensional vector.

Although the vectorized images are $p$ dimensional, the estimated intrinsic
dimension (ID) for the images is $4.18$ with a standard error of $0.1$,
determined using the two-nearest neighbor (NN) method (Facco et al., 2017).
This indicates that the high-dimensional vectorized images lie on a lower-
dimensional manifold, motivating the application of our proposed SkGP approach
to this data.

Out of $1667$ samples, we select every fourth sample point for the test set,
resulting in $n=1334$ training samples and $n_{new}=333$ test samples. As
depicted in Figure 9, the raw air quality monitor data displays
characteristics such as non-negativity, heavy-tailed distributions, non-
stationary patterns, and periodic behavior. To meet the normality assumption
for the error in (3), we apply a log transformation and standardize the
response, ensuring a mean of zero and a variance of one. The second row of
Figure 9 illustrates the log-transformed and standardized air quality data.
Many of the pixels in the satellite images are zero at all times and are
included to buffer the image to fit into a square. In our analysis, these
zeros are removed. The columns of the image predictor matrix, with zeros
removed, have dimension $1334\times 33068$ for the training data, and are
pixel-wise standardized to have zero mean and unit variance. We focus our
performance comparison on SkGP, BART and SkBART, considering them as the top
three competitors based on the simulation studies. While GP is also among the
top performers in the simulation studies, it is excluded from the comparison
due to its computational demands and memory intensity for this dataset. While
the data has a temporal component, our current analysis overlooks its time-
varying nature. Incorporating the temporal dynamics and capturing the evolving
associations between samples is a direction we plan to explore in future work.

Figure 9: Daily air quality monitor data in Las Vegas. Original (top) and
standardized for analysis (bottom). Data before 2019 contains irregular gaps
and is excluded. Figure 10: Near infrared image on July 2, 2019. The plot on
the left shows the original image. The right plot shows the same image with
screened out pixels in white. Interestingly, the independent screening
procedure selects contiguous chunks and borders in the image.

###  4.1 Results

The Nonparametric Independence Screening (NIS) method outlined in (Fan et al.,
2014) identifies $18640$ features out of $33068$ features which are marginally
related to the air quality. Figure 10 displays the pixels selected by NIS in a
representative multi-band image feature. Interestingly, even though the
screening procedure is independent for each pixel, contiguous patched of
pixels and boundaries around notable imaging patterns are screened out.

All competing models are implemented using $n=1334$ samples, where each sample
comprises air quality measurements as responses and $p=18640$ features.
Predictive inferences are generated for $n_{new}=333$ holdout samples. Figure
11 displays the point predictions and 95% predictive intervals for SkGP across
all time points, effectively capturing the trend in air quality responses.
Table 3 highlights the superior performance of SkGP compared to all
competitors, as evidenced by its lowest MSPE value. Although all competitors
exhibit under-coverage, potentially due to neglecting the time-varying nature
of the data, SkGP achieves the highest coverage (close to $80\%$) with
predictive intervals of comparable length to BART or SkBART. Overall, these
results underscore SkGP’s effectiveness in modeling the non-linear regression
relationship between air quality and multi-band satellite images.

Competitor | MSPE | Coverage | Length
---|---|---|---
SkGP | 0.327 | 0.784 | 1.165
BART | 0.369 | 0.739 | 1.300
SkBART | 0.536 | 0.613 | 1.159
Table 3: Mean squared Prediction Error (MSPE), length and coverage of 95%
predictive intervals for the competing methods SkGP, BART and SkBART for air
pollution data. Figure 11: Point prediction and 95% predictive interval at all
test samples of air pollution data for SkGP.

##  5 Conclusion and Future Work

This article is the first to present a novel Bayesian approach for predictive
inference of outdoor air quality using high-resolution satellite images, when
these images lie on a low-dimensional noise-corrupted manifold. Our
methodology exploits two powerful ideas, data sketching and stacking, to
eliminate the necessity for computationally demanding manifold structure
estimation, providing accurate point predictions and predictive uncertainties.
The computation of the posterior predictive distribution does not rely on MCMC
sampling, and our framework is amenable to parallel implementation, resulting
in substantial reductions in computation and storage costs. Empirical findings
underscore the significantly improved point prediction and predictive
uncertainty of our approach compared to existing methods. Future research
directions will extend our framework to handle large sample sizes using
distributed Bayesian inference (Guhaniyogi et al., 2022, 2023), exploring
applications to non-Gaussian or multivariate outcomes, and simultaneous
estimation of the intrinsic dimensionality of the manifold alongside
predictive inference for the outcome.

## References

  * Ahfock et al. (2017) D. Ahfock, W. J. Astle, and S. Richardson.  Statistical properties of sketching algorithms.  _arXiv preprint arXiv:1706.03665_ , 2017.
  * Apte et al. (2017) J. S. Apte, K. P. Messier, S. Gani, M. Brauer, T. W. Kirchstetter, M. M. Lunden, J. D. Marshall, C. J. Portier, R. C. Vermeulen, and S. P. Hamburg.  High-resolution air pollution mapping with google street view cars: exploiting big data.  _Environmental science & technology_, 51(12):6999–7008, 2017.
  * Arias-Castro et al. (2017) E. Arias-Castro, G. Lerman, and T. Zhang.  Spectral clustering based on local pca.  _The Journal of Machine Learning Research_ , 18(1):253–309, 2017.
  * Baraniuk et al. (2008) R. Baraniuk, M. Davenport, R. DeVore, and M. Wakin.  A simple proof of the restricted isometry property for random matrices.  _Constructive approximation_ , 28:253–263, 2008.
  * Belkin and Niyogi (2003) M. Belkin and P. Niyogi.  Laplacian eigenmaps for dimensionality reduction and data representation.  _Neural computation_ , 15(6):1373–1396, 2003.
  * Bleich et al. (2014) J. Bleich, A. Kapelner, E. I. George, and S. T. Jensen.  Variable selection for bart: an application to gene regulation.  2014\.
  * Breiman (1996) L. Breiman.  Stacked regressions.  _Machine learning_ , 24:49–64, 1996.
  * Breiman (2001) L. Breiman.  Random forests.  _Machine learning_ , 45:5–32, 2001.
  * Chen et al. (2010) M. Chen, J. Silva, J. Paisley, C. Wang, D. Dunson, and L. Carin.  Compressive sensing on manifolds using a nonparametric mixture of factor analyzers: Algorithm and performance bounds.  _IEEE Transactions on Signal Processing_ , 58(12):6140–6155, 2010.
  * Chen et al. (2018) Z. Chen, J. Fan, and R. Li.  Error variance estimation in ultrahigh-dimensional additive models.  _Journal of the American Statistical Association_ , 113(521):315–327, 2018.
  * Chipman et al. (2010) H. A. Chipman, E. I. George, and R. E. McCulloch.  Bart: Bayesian additive regression trees.  2010\.
  * Coifman and Lafon (2006) R. R. Coifman and S. Lafon.  Diffusion maps.  _Applied and computational harmonic analysis_ , 21(1):5–30, 2006.
  * Denison et al. (1998) D. G. Denison, B. K. Mallick, and A. F. Smith.  A bayesian cart algorithm.  _Biometrika_ , 85(2):363–377, 1998.
  * Dinh and Ho (2020) V. C. Dinh and L. S. Ho.  Consistent feature selection for analytic deep neural networks.  _Advances in Neural Information Processing Systems_ , 33:2420–2431, 2020.
  * Dobriban and Liu (2018) E. Dobriban and S. Liu.  A new theory for sketching in linear regression.  _arXiv preprint arXiv:1810.06089_ , 2018.
  * Drineas et al. (2011) P. Drineas, M. W. Mahoney, S. Muthukrishnan, and T. Sarlós.  Faster least squares approximation.  _Numerische mathematik_ , 117(2):219–249, 2011.
  * Drineas et al. (2012) P. Drineas, M. Magdon-Ismail, M. W. Mahoney, and D. P. Woodruff.  Fast approximation of matrix coherence and statistical leverage.  _The Journal of Machine Learning Research_ , 13(1):3475–3506, 2012.
  * Environmental Protection Agency (1996) Environmental Protection Agency.  _Air Quality Criteria for Particulate Matter_.  EPA/600/P-95/001aF. Office of Research and Development, Washington DC, 1996.
  * Facco et al. (2017) E. Facco, M. d’Errico, A. Rodriguez, and A. Laio.  Estimating the intrinsic dimension of datasets by a minimal neighborhood information.  _Scientific reports_ , 7(1):12140, 2017.
  * Fan et al. (2014) J. Fan, Y. Ma, and W. Dai.  Nonparametric independence screening in sparse ultra-high-dimensional varying coefficient models.  _Journal of the American Statistical Association_ , 109(507):1270–1284, 2014.
  * Gramacy (2020) R. B. Gramacy.  _Surrogates: Gaussian process modeling, design, and optimization for the applied sciences_.  CRC press, 2020.
  * Guhaniyogi and Dunson (2015) R. Guhaniyogi and D. B. Dunson.  Bayesian compressed regression.  _Journal of the American Statistical Association_ , 110(512):1500–1514, 2015.
  * Guhaniyogi and Dunson (2016) R. Guhaniyogi and D. B. Dunson.  Compressed gaussian process for manifold regression.  _The Journal of Machine Learning Research_ , 17(1):2472–2497, 2016.
  * Guhaniyogi and Scheffler (2021) R. Guhaniyogi and A. Scheffler.  Sketching in bayesian high dimensional regression with big data using gaussian scale mixture priors.  _arXiv preprint arXiv:2105.04795_ , 2021.
  * Guhaniyogi et al. (2022) R. Guhaniyogi, C. Li, T. D. Savitsky, and S. Srivastava.  Distributed bayesian varying coefficient modeling using a gaussian process prior.  _The Journal of Machine Learning Research_ , 23(1):3642–3700, 2022.
  * Guhaniyogi et al. (2023) R. Guhaniyogi, C. Li, T. Savitsky, and S. Srivastava.  Distributed bayesian inference in massive spatial data.  _Statistical science_ , 38(2):262–284, 2023.
  * Halko et al. (2011) N. Halko, P.-G. Martinsson, and J. A. Tropp.  Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions.  _SIAM review_ , 53(2):217–288, 2011.
  * Heffernan et al. (2023) C. Heffernan, R. Peng, D. R. Gentner, K. Koehler, and A. Datta.  A dynamic spatial filtering approach to mitigate underestimation bias in field calibrated low-cost sensor air pollution data.  _The Annals of Applied Statistics_ , 17(4):3056–3087, 2023.
  * Huang (2018) Z. Huang.  Near optimal frequent directions for sketching dense and sparse matrices.  In _International Conference on Machine Learning_ , pages 2048–2057. PMLR, 2018.
  * Jensen et al. (2021) K. Jensen, T.-C. Kao, J. Stone, and G. Hennequin.  Scalable bayesian gpfa with automatic relevance determination and discrete noise models.  _Advances in Neural Information Processing Systems_ , 34:10613–10626, 2021.
  * Lawrence and Hyvärinen (2005) N. Lawrence and A. Hyvärinen.  Probabilistic non-linear principal component analysis with gaussian process latent variable models.  _Journal of machine learning research_ , 6(11), 2005.
  * Le and Clarke (2017) T. Le and B. Clarke.  A bayes interpretation of stacking for m-complete and m-open settings.  2017\.
  * LeBlanc and Tibshirani (1996) M. LeBlanc and R. Tibshirani.  Combining estimates in regression and classification.  _Journal of the American Statistical Association_ , 91(436):1641–1650, 1996.
  * Li et al. (2022) D. Li, M. Mukhopadhyay, and D. B. Dunson.  Efficient manifold approximation with spherelets.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 84(4):1129–1149, 2022.
  * Li and Gu (2017) G. Li and Y. Gu.  Restricted isometry property of gaussian random projection for finite set of subspaces.  _IEEE Transactions on Signal Processing_ , 66(7):1705–1720, 2017.
  * Liu et al. (2021) Y. Liu, V. Rockova, and Y. Wang.  Variable selection with abc bayesian forests.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 83(3):453–481, 2021.
  * Maggioni et al. (2016) M. Maggioni, S. Minsker, and N. Strawn.  Multiscale dictionary learning: non-asymptotic bounds and robustness.  _The Journal of Machine Learning Research_ , 17(1):43–93, 2016.
  * Mahoney et al. (2011) M. W. Mahoney et al.  Randomized algorithms for matrices and data.  _Foundations and Trends® in Machine Learning_ , 3(2):123–224, 2011.
  * Mukhopadhyay and Dunson (2020) M. Mukhopadhyay and D. B. Dunson.  Targeted random projection for prediction from high-dimensional features.  _Journal of the American Statistical Association_ , 115(532):1998–2010, 2020.
  * Paciorek et al. (2008) C. J. Paciorek, Y. Liu, H. Moreno-Macias, and S. Kondragunta.  Spatiotemporal associations between goes aerosol optical depth retrievals and ground-level pm2. 5.  _Environmental science & technology_, 42(15):5800–5806, 2008.
  * Pavlyshenko (2020) B. M. Pavlyshenko.  Using bayesian regression for stacking time series predictive models.  In _2020 IEEE Third International Conference on Data Stream Mining & Processing (DSMP)_, pages 305–309. IEEE, 2020.
  * Planet Team (2017) Planet Team.  Planet application program interface: In space for life on earth.  _San Francisco, CA_ , 2017(40):2, 2017.  URL https://www.planet.com.
  * Raftery et al. (1997) A. E. Raftery, D. Madigan, and J. A. Hoeting.  Bayesian model averaging for linear regression models.  _Journal of the American Statistical Association_ , 92(437):179–191, 1997.
  * Tenenbaum et al. (2000) J. B. Tenenbaum, V. d. Silva, and J. C. Langford.  A global geometric framework for nonlinear dimensionality reduction.  _science_ , 290(5500):2319–2323, 2000.
  * Thanei et al. (2017) G.-A. Thanei, C. Heinze, and N. Meinshausen.  Random projections for large-scale regression.  _Big and Complex Data Analysis: Methodologies and Applications_ , pages 51–68, 2017.
  * Titsias and Lawrence (2010) M. Titsias and N. D. Lawrence.  Bayesian gaussian process latent variable model.  In _Proceedings of the thirteenth international conference on artificial intelligence and statistics_ , pages 844–851. JMLR Workshop and Conference Proceedings, 2010.
  * Weingessel and Hornik (2000) A. Weingessel and K. Hornik.  Local pca algorithms.  _IEEE Transactions on neural Networks_ , 11(6):1242–1250, 2000.
  * Williams and Rasmussen (2006) C. K. Williams and C. E. Rasmussen.  _Gaussian processes for machine learning_ , volume 2.  MIT press Cambridge, MA, 2006.
  * Wolpert (1992) D. H. Wolpert.  Stacked generalization.  _Neural networks_ , 5(2):241–259, 1992.
  * Woodruff et al. (2014) D. P. Woodruff et al.  Sketching as a tool for numerical linear algebra.  _Foundations and Trends® in Theoretical Computer Science_ , 10(1–2):1–157, 2014.
  * Yang and Dunson (2016) Y. Yang and D. B. Dunson.  Bayesian manifold regression.  2016\.
  * Yao et al. (2018) Y. Yao, A. Vehtari, D. Simpson, and A. Gelman.  Using stacking to average bayesian predictive distributions (with discussion).  2018\.
  * Yao et al. (2022a) Y. Yao, G. Pirs, A. Vehtari, and A. Gelman.  Bayesian hierarchical stacking: Some models are (somewhere) useful.  _Bayesian Analysis_ , 17(4):1043–1071, 2022a.
  * Yao et al. (2022b) Y. Yao, A. Vehtari, and A. Gelman.  Stacking for non-mixing bayesian computations: The curse and blessing of multimodal posteriors.  _The Journal of Machine Learning Research_ , 23(1):3426–3471, 2022b.
  * Zhang (2004) H. Zhang.  Inconsistent estimation and asymptotically equal interpolations in model-based geostatistics.  _Journal of the American Statistical Association_ , 99(465):250–261, 2004.
  * Zhang et al. (2013) L. Zhang, M. Mahdavi, R. Jin, T. Yang, and S. Zhu.  Recovering the optimal solution by dual random projection.  In _Conference on Learning Theory_ , pages 135–157, 2013.
  * Zhao et al. (2018) J. Zhao, L. Chen, W. Pedrycz, and W. Wang.  Variational inference-based automatic relevance determination kernel for embedded feature selection of noisy industrial data.  _IEEE Transactions on Industrial Electronics_ , 66(1):416–428, 2018.