# SHIMM as an atmospheric profiler on the Nickel Telescope

Ollie Jackson  University of California Santa Cruz  Maaike A.M. van Kooten
National Research Council Canada Herzberg Astronomy and Astrophysics Research
Centre  Saavidra Perera  University of California San Diego  Rebecca Jensen-
Clem  University of California Santa Cruz  Phil Hinz  University of California
Santa Cruz

## Abstract

Optimal atmospheric conditions are beneficial for detecting exoplanets via
high contrast imaging (HCI), as speckles from adaptive optics’ (AO’s)
residuals can make it difficult to identify exoplanets. While AO systems
greatly improve our image quality, having access to real-time estimates of
atmospheric conditions could also help astronomers use their telescope time
more efficiently in the search for exoplanets as well as aid in the data
reduction process. The Shack-Hartmann Imaging Motion Monitor (SHIMM) is an
atmospheric profiler that utilizes a Shack-Hartmann wavefront sensor to create
spot images of a single star in order to reconstruct important atmospheric
parameters such as the Fried parameter ($r_{0}$), $C_{n}^{2}$ profile and
coherence time. Due to its simplicity, the SHIMM can be directly used on a
telescope to get in situ measurements while observing. We present our
implementation of the Nickel-SHIMM design for the one meter Nickel Telescope
at Lick Observatory. We utilize an HCIPy simulation of turbulence propagating
across a telescope aperture to verify the SHIMM data reduction pipeline as we
begin on-sky testing. We also used on-sky data from the AO system on the Shane
Telescope to further validate our analysis, finding that both our simulation
and data reduction pipeline are consistent with previously determined results
for the Fried parameter at the Lick Observatory. Finally, we present first
light results from commissioning of the Nickel-SHIMM.

Keywords:  astronomical instrumentation, atmospheric profiling, Fried
parameter, site testing, High Contrast Imaging

##  1 Introduction

Over 5000 exoplanets have been detected and confirmed since the first was
discovered in the early 1990’s, and exoplanet science was one of the most
largely discussed fields of research in the 2020 decadal survey [1]. From
exoplanets, we can learn about the requirements for planet formation, whether
the formation of our Solar System is unique, and we can pursue questions of
astrobiology and life among the cosmos. There are a number of methods of
exoplanet detection that allow astronomers to research these questions, one of
which is high contrast imaging (HCI). HCI involves directly taking images of
exoplanets through the use of a coronagraph that blocks the light from the
host star, allowing for dim objects such as exoplanets to be visible.

For ground-based HCI, the instrument’s performance is highly dependent on the
atmospheric conditions. In Figure 1 from Cantalloube et. al. (2020)[2], the
same star is depicted on nights with different conditions, the left from a
night with good conditions and the right from a night with poor conditions. On
the right, the smearing and spreading of light around the star is referred to
as a wind-driven halo, and has been correlated to the atmospheric conditions
of high-altitude turbulent layers. The wind-driven halo structure makes it
challenging to search the region close to the star for exoplanets, and
demonstrates the necessity for optimal atmospheric conditions when using HCI
to detect exoplanets with small orbital radii.

Figure 1:  A night of good atmospheric conditions on the left, and poor
atmospheric conditions on the right, with the wind-driven halo effect clear on
the right. Taken from Cantalloube [2]

Knowing that the atmospheric conditions are sub-optimal just before starting a
long exposure on a target would allow astronomers to use their time on science
objects that may require less than perfect conditions (for example, searching
for an exoplanet further away from a host star on a night of poor conditions,
or even focusing on a companion close to a star on a night of excellent
conditions). Understanding the on-site atmospheric conditions is then
essential for ground-based astronomy because it will allow astronomers to use
their telescope time more effectively, and do better science as a result. This
has applications beyond just HCI, and can also be applied to queue-based
observing such as that done by ESO at the VLT.

There are a number of seeing monitors and atmospheric profilers that are used
frequently in astronomy to determine seeing conditions, namely DIMM
(Differential Image Motion Monitor [3]), MASS (Multi Aperture Scintillation
Sensor [4]),SCIDAR (SCIntillation Detection And Ranging [5]), and SLODAR
(SLOpe Detection And Ranging [6]). DIMM utilizes a mask with two subapertures
to get two images of a star. It then calculates the variance of the
differential motion between the image locations using the centroids to get an
estimate of the Fried parameter ($r_{0}$) [3]. While DIMM is a seeing monitor
that gives us a quantified estimate of the seeing conditions, MASS is an
atmospheric profiler that estimates a six-layer atmospheric turbulence
profile. By dividing the telescope aperture into rings, MASS can count the
number of photons on each ring to determine scintillation patterns and rebuild
turbulent strengths at defined heights above the telescope [4]. SCIDAR, like
MASS, utilises scintillation information to provide a higher-resolution
atmospheric profile by observing. However, unlike DIMM and MASS, it observes
two stars as opposed to just one. SCIDAR utilizes the notion that light from
two stars with some detectable angular separation must pass through the same
turbulence, such that similarly aberrated wavefronts reach the camera from two
different stars. A turbulence profile is recreated based on the correlation of
how each wavefront is distorted[5].

SLODAR analyzes where the light from two stars forms images after passing
through a Shack-Hartmann Wavefront Sensor (SHWFS) [6]. Described in further
detail in subsequent sections, a SHWFS is comprised of a lenslet array and
detector. When light passes through the lenslet array, it focuses into
multiple “spot images”, as visualized in the central row of Figure 2; an image
forms for every lenslet that light passes through. Since SLODAR looks at the
light from two stars, it creates two spot arrays (Figure 2) from the light of
each star, and correlations drawn from these spot arrays allow us to rebuild a
turbulence profile [7].

Each of these previously existing methods of estimating and measuring
turbulence profiles and/or measures of the seeing conditions have their
benefits and disadvantages. While all provide invaluable information to
astronomers, DIMM cannot rebuild a turbulence profile, and both SCIDAR and
SLODAR are necessarily large and expensive in order to produce their detailed
turbulence profiles. The Shack Hartmann Imaging Motion Monitor (SHIMM)
utilizes the same method as SLODAR to acquire $r_{0}$ and develops on the
functionality of DIMM by using a SHWFS instead of a mask to create
subapertures [8], while simultaneously being a low-cost and portable device
that can be more accessibly used by astronomers. The SHIMM’s simplicity allows
it to be used directly on a telescope to get in situ measurements while
observing.

Here we present our implementation of the SHIMM instrument and optical design
onto the 1-meter Nickel telescope (hereafter referred to as Nickel-SHIMM; see
Section 2.1), as well as the accompanying High Contrast Imaging Python
(HCIPy)[9] simulation and SHIMM data reduction pipeline/algorithm in Section
3. We then present the verification of our simulation and data pipeline,
followed by instrument commission and first results in Section 4.

##  2 SHIMM: Shack-Hartmann Imaging Motion Monitor

SHIMM is a portable, low-cost instrument that utilizes a single star to
calculate and estimate three important parameters that deepen our
understanding of atmospheric conditions at the telescope; the Fried parameter
($r_{0}$), the coherence time ($\tau_{0}$), and a three-layer profile of
atmospheric turbulence as a function of altitude ($C_{n}^{2}$ profile) [8].
The SHIMM utilizes the same methodology as the SLODAR to estimate $r_{0}$ from
subapertures of a single star, while SLODAR and SHIMM use different
methodologies to calculate $\tau_{0}$ and the $C_{n}^{2}$ profile. SLODAR uses
the correlation of spots from two stars for the $C_{n}^{2}$ profile while
SHIMM uses only one, as visualized in Figure 2 (noting that SLODAR uses a
single star for its estimation of $r_{0}$).

Figure 2: A comparison of the SLODAR and SHIMM designs

Beyond just the physical instrument, SHIMM uses a model of the instrument and
a data pipeline that reconstructs the above parameters. In this paper, we
focus on SHIMM’s calculation of $r_{0}$ as a measure of atmospheric
conditions, and in future research we will describe SHIMM’s estimation of
$\tau_{0}$ and $C_{n}^{2}$ profile. The $r_{0}$ value is the length over which
the phase of the wavefront changes by 1 radian. This parameter then gives a
quantitative description of how turbulent the atmosphere is, and serves as a
measure of the atmospheric conditions.

Since the SHIMM simply utilizes a SHWFS, it can be easily installed onto a
telescope (by redirecting some of the science light to the SHIMM) or by using
open loop adaptive optics (AO) wavefront sensing data. This provides a unique
opportunity to provide atmospheric information directly from the telescope
rather than being installed nearby or on an adjacent mountain. As a result,
the information we get from SHIMM will directly reflect the behavior of the
atmosphere seen by the telescope and gives astronomers an excellent
understanding of the on-site atmospheric conditions. Below we outline our
implementation of SHIMM on the Nickel telescope (Nickel-SHIMM), which utilizes
the same concepts and algorithms as SHIMM to measure $r_{0}$, $\tau_{0}$, and
the $C_{n}^{2}$ profile.

###  2.1 The Nickel-SHIMM Instrument Design

Figure 3: Diagram of the optical design for the Nickel telescope, including
visuals of the equipment used. The Nickel’s optical pathway is shown in red,
and Nickel-SHIMM’s in yellow demonstrating functionality of the flip-mirror.

The Nickel served as an optimal telescope to demonstrate SHIMM on a larger,
observationally driven telescope due to its aperture size and the resulting
subaperture size of our lenses. The subaperture size needs to be roughly near
the calculated $r_{0}$ value, and anywhere between 5 cm and 15 cm is a
reasonable range for the on-site $r_{0}$. The Nickel’s 40-inch diameter (1.016
m) with 12 subapertures across our lenslet array means that the diameter of
each of our subapertures is around 8 cm (1.016 m / 12 subapertures), so it
falls well within the 5 cm to 15 cm range of reasonable $r_{0}$ values.
SHIMM’s original testing (described in Perera et. al.[8]) utilized telescope
apertures of 9.25 and 11 inches, but the larger the telescope aperture, the
more subapertures are available, allowing for a more accurate calculation of
$r_{0}$ with more values to average together[8]. An even larger telescope than
the Nickel would have provided us with more subapertures, but the Nickel also
serves as an ideal testing site because it’s slightly under-subscribed. The
flexibility in the Nickel observing schedule allowed for us to propose
engineering time that would have minimal impact on other science being done.

The Nickel’s proximity to the Shane telescope further allows us to compare
results from the Nickel to the Shane, and test Nickel-SHIMM on calibration
data from Shane’s AO system (see Section 3.4). This will serve as an excellent
check of our results, and is another reason why we proposed to implement SHIMM
on the Nickel telescope. In addition, the Lick Observatory is easily
accessible from our research institute, allowing for easy travel to the site
for installation or in the event of a malfunction.

Our physical implementation of the Nickel-SHIMM instrument on the 1-m Nickel
Telescope consists of three main elements; the flip mirror, the Shack-Hartmann
lenslet array, and the camera. The flip mirror is an off-the-shelf ThorLabs
MFF102, which can be remotely rotated $90^{\circ}$ in or out of the optical
pathway to direct light out of the Nickel’s pathway and into our camera. For
our instrument, we used the 12.5 mm (0.5 inch) diameter mirror due to limited
space within the TUB. The lenslet array (MLA300-14AR from ThorLabs) contains
square lenses with a diameter of 295 $\mu$m (pitch of 300 $\mu$m) and focal
lengths of 14.3 mm, and is placed in the shaft of our CMOS sensor camera (a
Mako G-040 GigE camera from Allied Vision). Figure 3 displays our equipment
and a general diagram of how our instrument is set up.

The camera and flip mirror are operated via an Intel NUC computer (mounted
onto the telescope). To facilitate the power over ethernet required by the
camera, we have a smart network switch able to supply power to the camera
since it is unsupported by the NUC. This allows us to remotely turn on and off
the camera as well, making remote observation possible.

As mentioned briefly in Section 1, SHIMM utilizes a SHWFS, which consists of a
lenslet array and detector (camera). A lenslet array produces multiple images
from which we calculate $r_{0}$. Each lens in the lenslet array samples a
small portion of the incoming wavefront and creates its own spot in the focal
plane which is called a subaperture. A flat, undistorted (light) wave passing
through a lenslet array would form subaperture images on the focal plane at
the center of each lenslet, but when the incoming lightwaves are distorted by
atmospheric turbulence the images form off-center. The camera (detector) lies
in the focal plane of the lenslet array and captures where the subaperture
images form. By examining how a subaperture spot moves we can get the local
gradient value for that part of the aberrated wavefront. The phase shift of
the aberrated wavefront (i.e. how much the incoming lightwave has been
distorted from a plane wave) can be reconstructed from these slopes/gradients.

In addition to the flip mirror and SHWFS, our instrument contains a
collimating lens (focal length of 60mm) between the flip mirror and lenslet
array, which ensures that the lightwaves entering the SHWFS are parallel and
reimages the pupil plane for the SHWFS. There is also a transfer optic lens
placed before the flip mirror, which is part of the existing Nickel guider
optics, and which we then had to account for and incorporate into our design.
The transfer optic prescription is undocumented, but we estimate it to have a
focal length between 170 and 200 mm based on on-site testing. There is some
uncertainty in the focal length of the transfer optic, but it only changes the
final location of our pupil plane on the order of a millimeter. We find the
best pupil plane during alignment by evenly illuminating the telescope
aperture (via the dome flat lamps) to find the sharpest image of the telescope
aperture on our camera.

###  2.2 SHIMM on the Nickel

Between the nights of August 14 and August 15, 2022 the instrument was
installed onto the 40-inch Nickel Telescope located at the Lick Observatory on
Mount Hamilton. Our flip mirror, collimating lens, and camera (containing our
Shack-Hartmann lenslet array) were installed in the TUB of the Nickel
telescope (which houses the guider arm) at the base of the Cassagrian design,
before the Nickel’s own CCD camera (Figure 4)111For the Nickel Telescope’s
hardware overview, see
mthamilton.ucolick.org/techdocs/telescopes/Nickel/hw_overview/. We adjusted
the location of the camera using Vega as a bright star guide to ensure it was
centered on the optical path and images from the camera were in focus.
Additionally, our NUC computer containing the software necessary to operate
the camera and flip mirror was mounted to the outside of the telescope.

Figure 4: Physical Nickel-SHIMM instrument installed in the TUB of the Nickel
telescope. Yellow line indicates Nickel-SHIMM’s optical path.

We operate the camera via a Python interface, saving data cubes containing
between 500 and 2000 frames and using short exposures of 1-2 ms. The spots
images attained are consistent with our simulated predictions, and in a series
of frames the effects of turbulence can clearly be seen in the spots (both
through movement of the spot centroids and fluctuations in the brightness).
Following the procedure outlined in Section 3, the packets of images are
easily passed through the SHIMM algorithm to calculate $r_{0}$.

##  3 HCIPy and SHIMM Algorithm

###  3.1 HCIPy Simulation

We started with an end-to-end simulation of our optical system using the High
Contrast Imaging Python (HCIPy) [9] package to simulate light propagating
through changing turbulence before passing through the telescope. HCIPy allows
us to incorporate the specifications of the telescope (i.e. the diameter of
the primary mirror, the central obstruction due to the secondary mirror, the
science wavelength of light, and the throughput). For our simulation we
assumed a stellar magnitude of 4 and a throughput of 0.5. We additionally
assumed an outer scale ($L_{0}$) of 50 meters and a wind speed of 6 meters per
second in both the x and y directions (overall wind speed of roughly 8.5
meters per second). Based on the quantum efficiency of our camera, we chose a
wavelength of 500 nm 222See instrument specifications from Allied Vision at
www.alliedvision.com/en/camera-selector/detail/mako/g-040/. We incorporate our
SHWFS into the simulation, including the number of lenslets and subaperture
size, the focal length of the lenslets, and the lenslet array pitch (see
Section 2.1). By including the SHWFS in our simulation, we are able to
simulate the same subaperture spot images we will get on sky. We define a
single frozen-flow atmospheric layer from a set $r_{0}$ using HCIPy’s
implementation of an infinite von Karman layer. Finally, we propagate the
wavefront through the defined atmospheric layer, and at each timestep capture
the image of the wavefront through the simulated SHWFS while evolving the
atmospheric turbulence. The result is a stack of simulated spot images, an
example of which is shown in the left of Figure 5.

Figure 5: From left to right; a frame from HCIPy simulation, template for spot
locations, and mask.

We must ensure that we remove any bias coming from the simulation setup
itself. To do this, we use HCIPy to model SHWFS images that haven’t been
distorted by turbulence, and subtract the resulting slopes off of all
simulated images. We start with a singular reference image of our wavefront
and a mask, which encodes the spot locations for an undistorted wavefront and
maximum spot displacement. We also create a theoretical template of where the
images would form from an undistorted wavefront passing through the lenslets.
The difference between the mask and the template is that the mask takes each
spot where an undistorted image forms and encodes the maximum displacement
from that spot based on the diameter of each subaperture (visualized in Figure
5). From these (the reference image, the mask, and the template) we locate the
spots and centroids of the reference image, and get bias slopes that we must
subtract off of our simulated data.

Getting the theoretical/reference slopes is a very similar process to getting
the bias slopes; we repeat the process with the wavefront propagating through
atmospheric turbulence as opposed to a singular reference wavefront. Our
theoretical slopes come from simulated data where we set $r_{0}$ to be the
subaperature diameter of the lenslets (the diameter of the telescope divided
by the number of lenslets across our lenslet array). We advance a wavefront
through atmospheric layers and take images at each timestep. We do this
multiple times with our atmosphere evolving to sufficiently estimate the
statics of the system.For every timestep, we get slopes that describe the
differential of where the centroid of the image is relative to where it’s
theoretically supposed to be on our template. These reference slopes, minus
the bias slopes, are used to calculate $r_{0}$ for each real data set.

###  3.2 Extraction of Fried Parameter ($r_{0}$)

We estimate $r_{0}$ from the SHWFS subapertures using the methodology outlined
in Perera et. al. (2023)[8] and in Butterley et. al. (2006)[7]. From a time-
averaged cube of SHWFS images, we get the measured auto-covariance of the
centroids. We additionally create a cube of simulated SHWFS images in the
HCIPy simulation, and the auto-covariance of the simulated spot centroids
serves as our theoretical model. Theoretical/reference and on-sky auto-
covariance maps are created for both the x and y directions. A visual
comparison of the auto-covariance maps in the x direction for reference and
real data is shown in Figure 6. We then fit the theoretical auto-covariance
and on-sky auto-covariance maps using a one-dimensional strip through each map
as shown in Equation 1; we use the central strip, as the signal-to-noise ratio
(SNR) decreases as we move away from the center [8]. A first estimation of
$r_{0}$ for each direction is then found via Equation 1, where the
$\frac{-3}{5}$ comes from Kolmogorov’s power law, $m$ is the fitted parameter
(ratio of the measured auto-variance to the theoretical auto-covariance), $s$
is the subaperture diameter of the lenslets in the Shack-Hartmann lenslet
array, and $i$ indicates x or y direction. The $r_{0,i}$ values for x and y
are then averaged to get $r_{0}$.

$r_{0i}=m_{i}^{\frac{-3}{5}}*s$ (1)

We must also take into consideration the airmass above the telescope in order
to get $r_{0}$. We do this using the Pyephem library, which (given the date
and time of our observation, location of the telescope, and star name)
provides the zenith angle of our star. In our analysis, we noted that
Pyephem’s date and time was an hour behind Universal Standard time, and had to
adjust our calculation of the zenith angle as a result. The airmass is given
by 1/cosine(zenith angle), and the final $r_{0}$ value considering airmass is
given in Equation 2, where $r_{0}$ in Equation 2 is the average of the two
$r_{0,i}$ values calculated via Equation 1.

$r_{0,\text{air}}=r_{0}*\text{airmass}^{\frac{-3}{5}}$ (2)

Figure 6: Reference covariance matrix compared to covariance matrix from on-
sky data. The central strips are compared to calculate $r_{0}$.

###  3.3 Verifying Nickel-SHIMM and HCIPy Simulation

In order to verify the accuracy of our simulation, we created multiple
simulated datasets using various $r_{0}$ values. The difference between these
simulated datasets and our the simulated data we used to get our reference
slopes is determined entirely by $r_{0}$; an $r_{0}$ value equal to the
subaperture diameter of the lenslets defines our theoretical values. We
created simulated data sets with $r_{0}$ values ranging from 7 cm to 25 cm,
then for each data we recalculated $r_{0}$ from the spot images as described
in the previous section.

Figure 7 shows the theoretical $r_{0}$ values used to create the data plotted
against the $r_{0}$ values recalculated using the SHIMM algorithm. Our
simulation produces realistic results up to around 18 cm, and then it begins
to stray from where we would expect. For our purposes, however, an $r_{0}$
value within one to two centimeters of the actual value is sufficient to
accurately estimate the atmospheric conditions, and therefore we consider our
simulation to be working correctly.

Figure 7: Reconstructed r0 values from SHIMM algorithm (blue dots) compared
with the ”real” input r0 values (in red)

###  3.4 Testing on Shane AO Data

We also had access to calibration data for the AO system on the Shane
telescope also located on Mount Hamilton. The Shane telescope has an AO system
that allows it to correct for the effects of atmospheric turbulence in science
images. When the Shane’s AO system is run in “open loop” (the deformable
mirror is not correcting for turbulence), its SHWFS is measuring the full
atmospheric turbulence. A short amount of this open loop data is taken for
calibration purposes when Shane’s AO is used for observing, and occasionally
this data gets saved, allowing us to test our SHIMM algorithm and calculate
$r_{0}$. We adjusted our HCIPy simulation and SHIMM algorithm to fit the
specifications of the Shane telescope, which also involved creating new
reference images, templates, and masks.

Figure 8: Histogram of calculated $r_{0}$ values overlaid with Kernel density
function (estimation of the probability density using a non-parametric model).
In red is the mean $r_{0}$ determined by Walters[10], which is consistent to a
$\sigma$ range (in grey) with the mean $r_{0}$ on the Shane (shown in blue).

Using approximately 15 nights of data we created a histogram of the calculated
$r_{0}$ values, and compared this with the average $r_{0}$ on Mount Hamilton
as determined previously in Walters (1997)[10]. The calculated mean from our
nights of data using our SHIMM algorithm was 10.1 cm $\pm$ 5.7 cm, and the
mean determined in Walters (1997)[10] was 11.5 centimeters [10]. The value
from Walters falls within a $\sigma$ standard deviation range around our mean;
our calculated values are consistent with other determinations of $r_{0}$ at
the Lick Observatory (see Figure 8).

The SHIMM algorithm’s ability to recalculate $r_{0}$ from simulated data to a
reasonable degree, and verify the average $r_{0}$ value for the Lick
Observatory using on-sky data, reinforces that both the simulation and
algorithm are working as expected and are ready for on-sky testing.

##  4 On-Sky Testing and Results

We were granted frequent observation time on the Nickel from August 2022
through July 2023, typically collecting about an hour of data each night just
before sunrise. The results presented here are from data collected between May
and July of 2023. We focused our analysis on the star Deneb, and the binary
system Alberio was used to calibrate our system and determine the pixel scale
(which we found to be roughly 2 arcseconds per pixel). The exposure and gain
of our camera was adjusted nightly, but data was consistently collected using
exposure times between 1-2 milliseconds.

Figure 9: Left: Image captured by our camera during first light. Right: Image
from HCIPy simulation of the Nickel.

We calculated $r_{0}$ on data cubes of 500 frames at a time, using the
procedure outlined above. Preliminary results are displayed alongside our
$r_{0}$ estimations for the Shane telescope in Figure 10. Across our nights of
data, we found the average $r_{0}$ value to be 11.6 cm $\pm$ 5.6 cm, which is
consistent with the mean $r_{0}$ of 11.5 cm determined by Walters (1997)[10],
and is further consistent with the results from testing Nickel-SHIMM on Shane
AO calibration data.

In Figure 7 we saw that above $r_{0}$ values of 18 cm that our SHIMM pipeline
began to overestimate the actual $r_{0}$ value. We stated that for our
purposes, this behavior would not drastically skew our results. This is
corroborated by the fact that our mean $r_{0}$ values are well below 18 cm,
and that the conditions at the Nickel rarely reach $r_{0}$ values above 18 cm.

Figure 10: Preliminary results of Nickel-SHIMM compared with tests of SHIMM
algorithm on the Shane telescope. Histograms for each group of results
overlaid with corresponding Kernel Density functions, with the mean $r_{0}$
value calculated in Walters (1997)[10] again shown in red.

##  5 Discussion and Conclusions

Our results for $r_{0}$ at the Nickel are statistically consistent with both
our tests of the SHIMM algorithm on Shane open-loop adaptive optics data, and
with the mean $r_{0}$ for the Lick Observatory in Walters (1997) [10]. This
gives us confidence that our implementation of Nickel-SHIMM is working
correctly as we continue testing and optimizing our analysis. While the
preliminary results for the Shane telescope and Nickel telescope are
consistent, the $r_{0}$ values calculated from Shane data display different
behavior than those from the Nickel telescope. There are a handful of reasons
for this, first that we had much more on-sky data for the Nickel telescope
than for the Shane telescope. It’s also important to note that our testing of
the SHIMM algorithm on Shane open-loop data only provided a test for our data
pipeline and simulation; we did not have timestamps or the star target for the
Shane data so we were unable to take the airmass above the telescope into
account. Finally, with the Shane data, we have the deformable mirror for which
the quality of the flatness might affect our measurements.

As visualized in Figure 9, our camera and lens are slightly off-axis,
resulting in a tilt in our images compared to our simulation. When on-sky
images were compared with the templates, the spots in the upper and left sides
consistently overlapped well with their theoretical locations, but the spots
in the bottom right side were on average further from their theoretical
locations. This likely results in an under-estimation of $r_{0}$, where
atmospheric conditions are determined to be worse than they really are,
because a number of subapertures are always abnormally far from their
theoretical, undistorted locations. In future analysis, these subapertures
will be cut when calculating $r_{0}$.

As mentioned, Nickel-SHIMM is able to provide a three-layer $C_{n}^{2}$
profile and $\tau_{0}$, in addition to $r_{0}$. These capabilities are
currently in development for our implementation of Nickel-SHIMM. These results
indicate that Nickel-SHIMM is operating as expected and is able to provide
accurate representations of the atmospheric conditions using in-situ
measurements. As expected, Nickel-SHIMM will be able to provide real-time
information of atmospheric conditions to astronomers in future
implementations.

##  6 Acknowledgements

This research was made possible in part from the National Science Foundation
(award no. 1852393) and the Lamat Fellowship for Undergraduate Research at the
University of California Santa Cruz (UCSC). We would also like to thank Philip
Rice for his key support of the on-site implementation of SHIMM. This work was
done in collaboration with the Center for Adaptive Optics and the Lab for
Adaptive Optics at UCSC. Thank you to Rebecca Jensen-Clem and Eleanor Gates,
who provided us with the open-loop data from the Shane AO system. Thank you to
UCO and the UCO Machine Shops at UCSC, particularly Pavl Zachary and all the
staff at Mt Hamilton who operate the telescopes there.

The University of California Santa Cruz is located on the unceded land of the
Awaswas speaking Uypi tribe, and is currently stewarded under the hands of the
Amah Mutsun Tribal Band. Mount Hamilton is located on the lands of the
Muwekma, Ohlone, and Tamien/Tamyen Nations. This research would not be
possible without the lands on which it was completed, and those lands belong
to their historic and current Indigenous peoples.

## References

  * [1] National Academy of Sciences, Engineering, . M., “Pathways to discovery in astronomy and astrophysics for the 2020s,” The National Academies Press (2021).
  * [2] Cantalloube, F., Farley, O., Milli, J., Bharmal, N., Bradner, W., Correia, C., Dohlen, K., Henning, T., Osborn, J., Por, E., Valdez, M. S., and Vigan, A., “Wind-driven Halo in High Contrast Images,” Astronomy and Astrophysics 638 (2020).
  * [3] Sarazin, M. and Roddier, F., “The eso differential image motion monitor,” Astronomy and Astrophysics 227, 294–300 (1990).
  * [4] Kornilov, V., Tokovinin, A., Vozyakova, O., Zaitsev, A., Potanin, S., and Sarazin, M., “Mass: a monitor of the vertical turbulence distribution,” in [Adaptive Optics System Technologies II ], Proc. SPIE 4893 (2003).
  * [5] Vernin, J. and Roddier, F., “Experimental determination of two-dimensional spatiotemporal power spectra of stellar light scintillation evidence for a multilayer structure of the air turbulence in the upper troposphere,” Journal of the Optical Society of America 63, 270–273 (1973).
  * [6] Wilson, R., “Slodar: measuring optical turbulence altitude with a shack-hartmann wavefront sensor,” Monthly Notices of the Royal Astronomical Society 337, 103–108 (2002).
  * [7] Butterley, T., Wilson, R., and Sarazin, M., “Determination of the Profile of Atmospheric Optical Turbulence Strength from SLODAR Data,” Monthly Notices of the Royal Astronomical Society 369, 835–845 (2006).
  * [8] Perera, S., Wilson, R. W., Butterley, T., Osborn, J., Farley, O. J. D., and Laidlaw, D. J., “SHIMM: a versatile seeing monitor for astronomy,” Monthly Notices of the Royal Astronomical Society 520, 5475–5486 (02 2023).
  * [9] Por, E. H., Haffert, S. Y., Radhakrishnan, V. M., Doelman, D. S., Van Kooten, M., and Bos, S. P., “High Contrast Imaging for Python (HCIPy): an open-source adaptive optics and coronagraph simulator,” in [Adaptive Optics Systems VI ], Proc. SPIE 10703 (2018).
  * [10] Walters, D. and Bradford, L., “Measurements of $r_{0}$ and $\theta_{0}$: Two decades and 18 sites,” Applied Optics 36 (1997).


# Redistribution in Public Project Problems via Neural Networks

Guanhua Wang  The University of AdelaideAdelaideSouth AustraliaAustralia5005
<EMAIL_ADDRESS>,  Wuli Zuo  The University of
AdelaideAdelaideSouth AustraliaAustralia5005<EMAIL_ADDRESS>and  Mingyu Guo  The University of AdelaideAdelaideSouth
AustraliaAustralia5005<EMAIL_ADDRESS>
(2021)

###### Abstract.

Many important problems in multiagent systems involve resource allocations.
Self-interested agents may lie about their valuations if doing so increases
their own utilities. Therefore, it is necessary to design mechanisms
(collective decision-making rules) with desired properties and objectives. The
VCG redistribution mechanisms are efficient (the agents who value the
resources the most will be allocated), strategy-proof (the agents have no
incentives to lie about their valuations), and weakly budget-balanced (no
deficits).

We focus on the VCG redistribution mechanisms for the classic public project
problem, where a group of agents needs to decide whether or not to build a
non-excludable public project. We design mechanisms via neural networks with
two welfare-maximizing objectives: optimal in the worst case and optimal in
expectation. Previous studies showed two worst-case optimal mechanisms for 3
agents, but worst-case optimal mechanisms have not been identified for more
than 3 agents. For maximizing expected welfare, there are no existing results.

We use neural networks to design VCG redistribution mechanisms. Neural
networks have been used to design the redistribution mechanisms for multi-unit
auctions with unit demand. We show that for the public project problem, the
previously proposed neural networks, which led to optimal/near-optimal
mechanisms for multi-unit auctions with unit demand, perform abysmally for the
public project problem. We significantly improve the existing networks on
multiple fronts: We conduct a GAN network to generate worst-case type profiles
and feed prior distribution into loss function to provide quality gradients
for the optimal-in-expectation objective. We adopt dimension reduction to
handle a larger number of agents and we adopt supervised learning into the
best manual mechanism as initialization, then leave it into unsupervised
learning. For the worst case, we get better results than the existing manual
mechanisms, and for the optimal-in-expectation objective, our mechanisms’
performances are close to the theoretical optimal performance.

public project, VCG redistribution mechanism, neural networks, GAN, dimension
reduction

††journalyear: 2021††copyright: acmcopyright††conference: IEEE/WIC/ACM
International Conference on Web Intelligence; December 14–17, 2021; ESSENDON,
VIC, Australia††booktitle: IEEE/WIC/ACM International Conference on Web
Intelligence (WI-IAT ’21), December 14–17, 2021, ESSENDON, VIC,
Australia††price: 15.00††doi: 10.1145/3486622.3493922††isbn:
978-1-4503-9115-3/21/12††ccs: Theory of computation Algorithmic mechanism
design

##  1\. Introduction

###  1.1. VCG Redistribution Mechanisms

Many important problems in multiagent systems are related to resource
allocations. The problem of allocating one or more resources among a group of
competing agents can be solved through economic allocation mechanisms that
take the agents’ reported valuations for the resources as input, and produce
an allocation of the resources to the agents, as well as payments to be made
by the agents. As a central research branch in economics and game theory,
mechanism design concerns designing collective decision-making rules for
multiple agents, to achieve desirable objectives, such as maximizing the
social welfare, while each agent pursues her own utility. A mechanism is
efficient if the agents who value the resource the most will get it. A
mechanism is strategy-proof if the agents have the incentives to report their
valuations truthfully, which is to say, an agent’s utility is maximized when
reporting her true valuation, no matter how the other agents report.

The Vickrey-Clarke-Groves (VCG) mechanism is a celebrated efficient and
strategy-proof mechanism. Under the VCG mechanism, each agent $i$ reports her
private type $\theta_{i}$. The outcome that maximizes the agents’ total
valuations is chosen. Every agent is required to make a VCG payment
$t(\theta_{-i})$, which is determined by the other agents’ types. An agent’s
VCG payment is often described as how much this agent’s presence hurts the
other agents, in terms of the other agents’ total valuations. The total VCG
payment may be quite large, leading to decreased welfare for the agents. In
particular, in the context of the public project problem, where the goal is
often to maximize the social welfare (the agents’ total utility considering
payments), having large VCG payments are undesirable.

To address the welfare loss due to the VCG payments, Cavallo (Cavallo, 2006)
suggested that we first execute the VCG mechanism and then redistribute as
much of the payments back to the agents, without violating the efficiency and
strategy-proofness of the VCG mechanism, and in a weakly budget-balanced way.
This is referred to as the VCG redistribution mechanism. The amount that every
agent receives (or pays additionally) is called the redistribution payment. To
maintain efficiency and strategy-proofness of VCG, the redistribution payment
of an agent is required to be independent of her own valuation. To maintain
weakly budget-balance, the total amount redistributed should never exceed the
total VCG payment. The redistribution payment is characterized by a
redistribution function $h$, where $h(\theta_{-i})$ represents agent $i$’s
redistribution payment.

There have been many successes on designing redistribution mechanisms for
various multi-unit/combinatorial auction settings (Cavallo, 2006; de Clippel
et al., 2014; Faltings, 2005; Guo and Conitzer, 2009; Moulin, 2009; Gujar and
Narahari, 2011; Guo, 2011, 2012; Guo and Conitzer, 2014), including a long
list of optimal/near-optimal mechanisms. On the other hand, there hasn’t been
comparable success in solving for optimal redistribution mechanisms for the
public project problem, despite multiple attempts (Naroditskiy et al., 2012;
Guo, 2016; Guo and Shen, 2017; Guo, 2019). In terms of optimal results,
Naroditskiy et al. (Naroditskiy et al., 2012) solved for the worst-case
optimal mechanism for three agents. Unfortunately, the authors’ technique does
not generalize to more than three agents. Guo (Guo, 2019) proposed a mechanism
that is worst-case optimal when the number of agents approaches infinity, but
for small number of agents, the mechanism is not optimal. For maximizing
expected welfare, there are no existing results, because it is difficult for
traditional mathematical analysis(eg: mixed integer programming) to maximize
the expectation of welfare.

###  1.2. Designing VCG Redistribution Mechanisms via Neural Networks

A recent emerging topic in mechanism design is to bring tools such as neural
networks from machine learning to design mechanisms (Duetting et al., 2019;
Golowich et al., 2018; Manisha et al., 2018; Shen et al., 2019; Wang et al.,
2021). Duetting et al. (Duetting et al., 2019) proposed a neural network
approach for the automated design of optimal auctions. They model an auction
as a multi-layer neural network and frame optimal auction design as a
constrained learning problem which can be solved using standard machine
learning pipelines. The training and testing type profiles are generated based
on the prior distribution. The cost function involves the mechanism objective
and the penalty for property violation. Essentially, neural networks were used
as tools for functional optimisation. Shen et al. (Shen et al., 2019) proposed
a neural network based framework to automatically design revenue optimal
mechanisms. This framework consists of a seller’s network, which provides a
menu of options to the buyers, and a buyer’s network, which outputs an action
that maximizes her utility. Wang et al. (Wang et al., 2021) studied mechanism
design for the public project problem and proposed several technical
innovations that can be applied to mechanism design in general to improve the
performance of mechanism design via neural networks.

Figure 1. Neural network structure reported by Manisha et al. (Manisha et al.,
2018)

The work by Manisha et al. (Manisha et al., 2018) is the first and only
attempt to design VCG redistribution mechanisms using neural networks. They
focused on multi-unit auctions with unit demand and studied both worst-case
and optimal-in-expectation objectives. By randomly generating a large number
of bid profiles, they train a neural network to maximize the total
redistribution payment, while enforcing that the total redistribution should
not exceed the total VCG payment. They modelled the redistribution function as
a simple network outlined in Figure 1. It is a fully connected network with
one hidden layer using ReLU activation. It takes the valuations of all the
agents other than agent $i$ herself as input, and outputs the predicted
redistribution payment for agent $i$. Their data is sampled randomly from
uniform distribution ($\theta_{i}\in Uniform(0,1)$).

###  1.3. Improved Neural Networks for Designing VCG Redistribution Mechanisms
for the Public Project Problem

In this paper, we train neural networks to design VCG redistribution functions
for the public project problem, which turns out to be a more challenging
setting compared to multi-unit auctions studied by Manisha et al. (Manisha et
al., 2018). The public project problem is a classic mechanism design problem
that has been studied extensively in economics and computer science (Mas-
Colell et al., 1995; Moore, 2006; Moulin, 1988). In this problem, $n$ agents
decide whether or not to build a non-excludable public project, for example, a
public bridge that can be accessed by everyone once built. Without loss of
generality, we assume that the cost of the project is $1$, and
$\theta_{i}(0\leq\theta_{i}\leq 1)$ is agent $i$’s valuation for the project
if it is built. If the decision is not to build, every agent retains her share
of the cost, which is $1/n$.

We first evaluate the simple multilayer perceptron (MLP) model proposed by
Manisha et al. (Manisha et al., 2018). That is, for each agent $i$, we train a
neural network that maps $\theta_{-i}$ to agent $i$’s redistribution. The
training and testing samples are randomly generated based on the prior
distribution. The cost function maximizes the mechanism design objective, as
well as enforces mechanism design constraints via penalty. We find that such a
simple MLP is not effective enough for the public project problem for the
following reasons:

  1. (1)

From our experiments, by randomly generating the type profiles, we are not
getting the true worst-case type profiles for the public project problem (it
is a coincidence that for multi-unit auctions with unit demand, it is a lot
easier to hit a worst case).

  2. (2)

Another challenge is the high input dimension when the number of agents is
large. For 100 agents, the neural network has to take a 99-dimensional input,
which is computationally unrealistic.

  3. (3)

In the public project problem, the agents’ collective payments differ
significantly between cases where the decision is to build the public project
and cases where the decision is not to build. The number of “build” samples in
a training batch significantly impacts the parameter gradient, which results
in a wild loss fluctuation during the training process.

To solve the aforementioned problems, we propose a novel neural network
approach to design redistribution mechanisms for the public project problems.
Our approach involves the following technical innovations.

##### GAN Network

For the worst-case objective, we introduce a generative adversarial network
(GAN) to generate worst-case type profiles, and then use these type profiles
to train the redistribution function. Our experiment shows that a mechanism
trained only using randomly generated data fails when facing type profiles
generated by GAN, so GAN is necessary and effective to derive the worst case.

##### Dimension Reduction

Instead of feeding $\theta_{-i}$ as input to train the mechanism, which has
$n-1$ dimensions, we extract a few expressive features from $\theta_{-i}$
(e.g. the maximum of types $\theta_{-i}$, the sum of $\theta_{-i}$ excluding
the maximum, etc.). This reduces the input dimension to 3. This helps the
neural network loss converge faster and still retain good performance.

##### Supervised Learning

Wang et al. (Wang et al., 2021) suggested that supervision to manual
mechanisms often outperforms random initialization in terms of training speed
by pushing the performance to a state that is already somewhat close to
optimality. In addition, unlike many other deep learning problems, for
mechanism design, there often exist simple and well-performing mechanisms that
can be used as starting points. In this particular problem, we first conduct
supervised learning to let the network mimic the state of art manual mechanism
(Guo, 2019), and then leave it to gradient descent. This approach saves time
for larger $n$ in our experiments.

##### Feeding Prior Distribution into Loss Function

We use probability density function (PDF) of the prior distribution to provide
quality gradients. In specific, for each valuation profile
$\theta=\\{\theta_{i}\\}(i=1..n)$ generated from a distribution $D$, we
randomly choose $\theta_{i}$ to be replaced by a randomly generated
$\theta^{\prime}_{i}$ from $Uniform(0,1)$ and update $\theta$ to be
$\theta^{\prime}$. This sample is then assigned a weight based on the PDF. In
experiments, we see that this approach significantly reduces the loss
fluctuation during training. One potential explanation (or observation) is
that this approach reduces the fluctuation in the proportion of “build” cases
among a batch.

With a more sophisticated network architecture due to the above technical
adjustments, we get better results for the worst-case than the state of the
art (Guo, 2019). For the optimal-in-expectation objective, our results are
close to the theoretical optimal values.

###  1.4. Organization

The rest of the paper is organized as followings: Section 2 gives a formal
description of the VCG redistribution mechanism for the public project
problem. Section 3 and 4 present our neural network approaches for the worst-
case and expectation optimal objectives respectively, with the explanation of
the network architecture, the technical adjustments, and loss function design.
Section 5 discusses the training of neural networks and the experimental
analysis. Section 6 concludes.

##  2\. Formal Model Description

For the public project problem, VCG redistribution mechanisms have the
following form (Naroditskiy et al., 2012):

  * •

Build the public project if and only if $\sum_{i}\theta_{i}\geq 1$.

  * •

If the decision is to build, agent i receives $\sum_{j\neq
i}\theta_{j}-h(\theta_{-i})$.

  * •

If the decision is not to build, then agent i receives
$(n-1)/n-h(\theta_{-i})$.

  * •

$h$ is an arbitrary function and $\theta_{-i}$ refers to the types from the
agents other than $i$ herself.

A VCG redistribution mechanism is characterized by the function $h$.

Due to Guo (Guo, 2019),

  * •

$S(\theta)=max\\{\sum_{i}\theta_{i},1\\}$ is exactly the first-best total
utility. (I.e., if the sum of types is higher than $1$, then the efficient
decision is to build. Otherwise, the efficient decision is not to build.)

  * •

The agents’ welfare (total utility considering payments) under type profile
$\theta$ is $nS(\theta)-\sum_{i}h(\theta_{-i})$, which is obtained via simple
algebraic simplification based on the definition of VCG redistribution
mechanisms.

We consider two objectives. One is to find a mechanism that maximizes the
worst-case efficiency ratio, and the other is to maximize the expected
efficiency ratio.

###  2.1. Worst-case Optimal Mechanism

The efficiency ratio $r$ is defined as the ratio between the achieved total
utility and the first best total utility:

$r=\frac{nS(\theta)-\sum_{i}{h(\theta_{-i})}}{S(\theta)}=n-\frac{\sum_{i}{h(\theta_{-i})}}{S(\theta)}$

The worst-case efficiency ratio is the worst case ratio between the achieved
total utility and the first best total utility, namely, the minimum of $r$
over all type profiles. Due to Guo (Guo, 2019), the mechanism has a worst-case
efficiency ratio $\alpha$ if and only if:

(1) $\forall\theta,(n-1)\leq\sum_{i}h(\theta_{-i})/S(\theta)\leq(n-\alpha)$

In Inequality 1, the left side is the constraint for weakly budget-balance,
and the right side corresponds to the definition of $\alpha$.

Therefore, taking the worst-case ratio as the objective, we need to design an
$h$ function that:

[3] α ∀θ,(n-1) ≤∑_ih(θ_-i)/S(θ)≤(n-α)

###  2.2. Optimal-in-Expectation Mechanism

For this objective, we maximize the expected efficiency ratio $r$ to 1, which
is equivalent to minimize $\sum_{i}{h(\theta_{-i})}/S(\theta)$ from above to
$n-1$, with the consideration of the weakly budget-balance constraint.

We are designing an $h$ function that:

[3] ¯∑_ih(θ_-i)/S(θ) ∀θ,(n-1)≤∑_ih(θ_-i)/S(θ)

##  3\. Worst-case Optimal Mechanism

In this section, we focus on the worst-case optimal mechanism. We first
describe our neural network approach by explaining the network architecture
and details of the relating techniques. Then we define the loss function.

###  3.1. Network Architecture

Figure 2. Worst-case optimal network architecture

We construct a neural network to determine the $h$ function. As illustrated in
Figure 2, it is a network system in which a GAN and an MLP interacting with
each other, we call it GAN+MLP.

The GAN works as a Generative Model and is used to generate special samples
(type profiles). It takes $n$ randomly generated values as input ideas, and
the output is $\theta=\\{\theta_{i}\\}(i=1...n)$. It is a fully connected
network with 4 hidden layers, and each hidden layer contains 64 nodes. For a
given batch ($batch\\_size=b$), the GAN generates $b$ type profiles:
$\theta^{(j)}\in batch=\\{\theta^{(1)},\theta^{(2)},...,\theta^{(b)}\\}$, with
the aim to maximize the difference between the maximum and the minimum of
$\sum_{i}{h(\theta^{(j)}_{-i})}/S(\theta^{(j)})$. It means to:

$maximize\
(\sum_{i}{h(\theta^{(j_{1})}_{-i})}/S(\theta^{(j_{1})})-\sum_{i}{h(\theta^{(j_{2})}_{-i})}/S(\theta^{(j_{2})}))$

where $\theta^{(j1)},\theta^{(j2)}\in batch$ is the sample that gives the
maximum and minimum of $\sum_{i}{h(\theta^{(j)}_{-i})}/S(\theta^{(j)})$,
respectively.

The MLP works as a Discriminative Model that learns the samples generated by
the GAN. The MLP is a fully connected neural network. For each agent $i$, the
network takes $\theta_{-i}$ as the input, and outputs the value of
$h(\theta_{-i})$. In the MLP, there are 6 hidden layers, each of which
contains 100 nodes and with ReLU as the activation function. We first train
the MLP under supervision to the best-performing manual mechanism and then
leave it to unsupervised learning. For unsupervised learning, our cost
function is the combination of design objective and also penalty due to
constraint violation.

###  3.2. Details of the Networks and Evaluations

To improve the result for the neural network, we adopt some technical
adjustments. We use a GAN instead of uniform to generate special cases in
order to find out the worst case. For cases with a greater number of agents
($n\geq 5$), we adopt two technical tricks: Dimension Reduction and Supervised
Learning.

####  3.2.1. GAN Network

In previous studies, the authors used random generation or fixed data to find
the worst case (Manisha et al., 2018).

As mentioned in Section 3.1, We propose a new GAN approach to find out the
worst case.

We conduct a contrast experiment to verify the validation of the GAN. We use
only data generated from uniform distribution to train a network, and test
this network with two sets of data.

  * •

Test Set A: 20000 data drawn from $Uniform(0,1)$

  * •

Test Set B: 10000 data generated from a trained GAN network and 10000 data
drawn from $Uniform(0,1)$

Figure 3 outlines the experimental results for $n=10$ showing the difference
of the the network performance on Test Set A and B. In the left figure, the
network is tested by randomly generated data set A. It gives $\alpha=0.896$,
and $\sum_{i}{h(\theta_{-i})}/S(\theta)$ is between 9 and 9.104. However, in
the right figure, the network performs poorly with significant violation of
the weakly budget-balance constraint ($\sum_{i}{h(\theta_{-i})}/S(\theta)$ is
from 3 to 9.2).

Figure 3. Spread of $\sum_{i}{h(\theta_{-i})}/S(\theta)$ evaluated using
random type profiles and GAN generated type profiles. Random generation fails
to generate true worst cases.

Therefore, random type profile generation as used in Manisha et al. (Manisha
et al., 2018) does not work for this problem. To get the worst-case
performance, we need a GAN network to generate higher quality worst-case
profiles, and then let the MLP learn these profiles.

####  3.2.2. Dimension Reduction

The MLP takes a $(n-1)$-dimensional input for $n$ agents, resulting in
expensive computation when $n$ is large. This motivates us to look for an
effective dimension-reducing technique.

We first manually design a list of features that describe $\theta_{-i}$, and
then experimentally search for a good combination of three features to be used
for dimension reduction purposes. (We essentially reduce $\theta_{-i}$ to
three dimensions this way.)

The features we consider include:

  * •

The highest type(s) from $\theta_{-i}$

  * •

The the lowest type(s) from $\theta_{-i}$

  * •

The sum of some types from $\theta_{-i}$

  * •

The standard deviation of some types from $\theta_{-i}$

  * •

The largest jump of adjacent types from $\theta_{-i}$
Here $jump\ of\ adjacent\ types$ is defined as: for a sorted list
${\theta_{-i}}$, there is $jump_{j}$ between $\theta_{j}$ and $\theta_{j+1}$:

$jump_{j}=\theta_{j+1}-\theta_{j}$

We first experimentally derive that the following features are more important
than the rest (i.e., removing them results in significant performance loss):

  * •

The highest type from $\theta_{-i}$

  * •

The the lowest type from $\theta_{-i}$

  * •

The sum of some type from $\theta_{-i}$

We then experimentally evaluate different combinations of the above features:

  1. (1)

Combination 1: the highest type(s) from $\theta_{-i}$ & the sum of all the
other types

  2. (2)

Combination 2: the highest type(s) from $\theta_{-i}$ & the difference between
the highest and the lowest type from $\theta_{-i}$

  3. (3)

Combination 3: the highest type(s) from $\theta_{-i}$ & the standard deviation
of all types from $\theta_{-i}$

  4. (4)

Combination 4: the highest type(s) from $\theta_{-i}$ & the standard deviation
of all the other types

  5. (5)

Combination 5: the highest type(s) from $\theta_{-i}$ & the largest jump of
adjacent types

  6. (6)

Combination 6: the highest type(s) from $\theta_{-i}$, the lowest type from
$\theta_{-i}$, & the largest jump of adjacent types

  7. (7)

Combination 7: the highest type(s) from $\theta_{-i}$, the sum of all the
other types & the largest jump of adjacent types

  8. (8)

Combination 8: the highest type(s) from $\theta_{-i}$, the lowest type from
$\theta_{-i}$, & the sum of all the other types

Figure 4 shows $\alpha$ of the networks trained with the different input
combinations against the number of agents $n$. It is found that Combination 1,
7 and 8 performs better than the other combinations.

Figure 4. The worst-case ratio $\alpha$ of the models trained with the input
generated by using different feature-combinations against the number of agents
$n=5,...,10$.

The worst-case ratio $\alpha$ of the models trained with the input generated
by using different feature-combinations against the number of agents
$n=5,...,10$.

The above dimension-reducing mechanism improves both the training speed and
sometimes the performance. We can infer that by using this dimension-reducing
mechanism, the training speed would have a more significant improvement with
the increase of the agent number $n$.

Figure 5. The loss during the training when using (left) and not using (right)
the dimension-reducing mechanism for $n=10\ (\theta_{i}\in Uniform(0,1))$.
Dimension reduction leads to faster convergence.

Figure 5 shows the difference of the loss of the network between using and not
using dimension-reducing mechanism for $n=10$. In the left figure, the
training and test loss is stabilized within 1000 iterations with the
application of the dimension-reducing mechanism, while in the right figure,
the losses still vibrate in a wider range till 2000 iterations.

Figure 6. Spreading of $\sum_{i}h(\theta_{-i})/S(\theta)$ when using (left)
and not using (right) the dimension-reducing mechanism for $n=10\
(\theta_{i}\in Uniform(0,1))$, data size = 10000, bins = 500. Dimension
reduction leads to better-performing mechanism.

Figure 6 shows the difference on the spreading of
$\sum_{i}{h(\theta_{-i})}/S(\theta)$ between using and not using dimension-
reducing mechanism for $n=10$. In the left figure, when using the dimension-
reducing mechanism, we get an expected performance of
$\sum_{i}{h(\theta_{-i})}/S(\theta)=9.003$, which is very close to $n-1=9$.
This indicates the corresponding mechanism is close to being optimal. (We
recall that for the expected performance, we want the expected ratio to be as
close to $n-1$ as possible.) When not using the dimension-reducing mechanism,
our performance is $9.377$, which is further away from $n-1=9$.

####  3.2.3. Supervised Learning

For a greater number of agents, it takes gradient descent a long time to fix
the constraint violations. Wang et al. (Wang et al., 2021) suggested that we
first supervise the neural network into the existing best manual mechanism,
and then leave it to unsupervised learning. With the best manual mechanism as
the starting point, better mechanisms can usually be found. The existing best
manual mechanism for worst-case optimal objective is reported by Guo (Guo,
2019).

###  3.3. Loss Function

####  3.3.1. GAN Network Loss

The GAN network has the following loss:

$\displaystyle loss_{GAN}=$ $\displaystyle
min\\{\sum_{i}{h(\theta^{(j_{1})}_{-i})}/S(\theta^{(j_{1})})\\}$
$\displaystyle-
max\\{\sum_{i}{h(\theta^{(j_{2})}_{-i})}/S(\theta^{(j_{2})})\\},\ $
$\displaystyle\theta^{(j_{1})},\theta^{(j_{2})}\in batch$

####  3.3.2. Supervised Loss

In supervised learning, we want the predicted $h$ to be as close as possible
to the best manual value $h\\_manual$ (Guo, 2019), so the loss is:

$loss_{supervised}=(h-h\\_manual)^{2}$

####  3.3.3. Unsupervised Loss

In the unsupervised learning stage, we want all the
$\sum_{i}{h(\theta_{-i})}/S(\theta)$ to maintain the weakly budget-balance
constraint, and let $\alpha$ to be close to 1, which means to make the upper
bound of $\sum_{i}{h(\theta_{-i})}/S(\theta)$ as small as possible.

The loss function includes two parts:

  * •

$objective\\_loss$ =
$(relu(\sum_{i}{h(\theta_{-i})}-(n-up\\_bound)S(\theta)))^{2}$

  * •

$constraint\\_loss$ = $(relu((n-1)S(\theta)-\sum_{i}{h(\theta_{-i})}))^{2}$

Since the weakly budget-balance is a strict constraint, while the objective is
soft, we add a multiplier $\epsilon$ to weaken the effect of
$objective\\_loss$. For the worst-case optimal network, we get the best
$\epsilon=0.01$ through experiments.

$\displaystyle loss_{unsupervised}$ $\displaystyle=\epsilon\cdot
objective\\_loss+constraint\\_loss$
$\displaystyle=objective\\_loss/100+constraint\\_loss$

##  4\. Optimal-in-Expectation Mechanism

We design the optimal-in-expectation mechanism with slight modifications based
on the worst-case mechanism.

The architecture of the MLP stays the same. We do not need a GAN for
generating the worst-case since the worst-case does not matter. For a large
agent number $n$, we also adopt dimension reduction and supervised learning as
we do for the worst-case mechanism. In addition, we feed the prior
distribution into the loss function to achieve a high-quality gradient. This
network is called MLP+FEED.

###  4.1. Feed Prior Distribution into Loss Function

In this problem, the decision to build or not to build significantly affects
the expectation of a training batch. For batches with different amounts of
“build” cases, the gradients fluctuate significantly, causing worse training
results and speed.

Wang et al. (Wang et al., 2021) discovered a way to insert the cumulative
distribution function (CDF) of the prior distribution into the cost function
to get more accurate loss function. The approach was shown to be effective for
optimal-in-expecation mechanism design. We adopt a similar idea, but use
probability density function (PDF) from the prior distribution to provide
quality gradients for our training process. For each valuation profile
$\theta=\\{\theta_{i}\\}(i=1..n)$ generated from a distribution $D$, we
randomly choose $\theta_{i}$ to be replaced by $\theta^{\prime}_{i}$ which is
regenerated from $Uniform(0,1)$ and update $\theta$ to be $\theta^{\prime}$.

The probability of the profile $\theta^{\prime}$ is proportional to
$PDF_{D}(\theta^{\prime}_{i})$, so the loss should also be multiplied by
$PDF_{D}(\theta^{\prime}_{i}$).

$loss_{unsupervised\\_feeding}=loss_{unsupervised}\cdot
PDF_{D}(\theta^{\prime}_{i})$

Here,

$PDF_{D}(\theta^{\prime}_{i})=10^{log\\_prob(\theta^{\prime}_{i})}$

$log\\_prob$ is provided by PyTorch to calculate PDF. PyTorch distributions
package is based on Schulman (Schulman et al., 2015). Our experiments show the
difference between feeding and not feeding the distribution into loss
function. Figure 7 shows that for normal distribution ($Normal(0.5,0.1),\
n=3$), feeding distribution into the loss function helps get a better gradient
and thus dramatically improves the test result.

Figure 7. Test loss (the distance from $\sum_{i}{h(\theta_{-i})}/S(\theta)$ to
$(n-1)$) during the training for distribution: $Normal(0.5,0.1),\ n=3)$.
Feeding prior distribution leads to faster convergence and better mechanism.

###  4.2. Loss Function

The loss for the optimal-in-expectation network is similar to that of the
worst-case MLP.

####  4.2.1. Supervised Loss

$loss_{supervised}=(h-h\\_manual)^{2}$

####  4.2.2. Unsupervised Loss

We use the square loss to approximate the loss for the expectation, since the
derivative of square loss ($f(x)=x^{2}$) is a linear function, and can
represent the strength of the gradient decedent, which is also linear
($f^{\prime}(x)=ax+b$).

Loss function consists of $objective\\_loss$ and $constraint\\_loss$, where

  * •

$objective\\_loss=(relu(\sum_{i}{h(\theta_{-i})}-(n-1)S(\theta)))^{2}$

  * •

$constraint\\_loss=(relu((n-1)S(\theta)-\sum_{i}{h(\theta_{-i})}))^{2}$

$objective\\_loss$ aims to push $\sum_{i}{h(\theta_{-i})}/S(\theta)$ close to
$(n-1)$, so that the mechanism redistributes as much of the collected VCG
payment back to the agents. $constraint\\_loss$ is for weakly budget-balance,
to make the total redistributed amount less than the total VCG payment.

Similar with the worst-case optimal network, we use a multiplier $\epsilon$ to
soften $objective\\_loss$, and the experimentally best $\epsilon=10^{-4}$. So,

$\displaystyle loss_{unsupervised}$ $\displaystyle=\epsilon\cdot
objective\\_loss+constraint\\_loss$
$\displaystyle=objective\\_loss/10000+constraint\\_loss$

With the feeding of prior distribution to loss function as described in
Section 4.1, we have:

$\displaystyle loss_{unsupervised\\_feeding}$
$\displaystyle=(objective\\_loss/10000+constraint\\_loss)*PDF_{D}(\theta^{\prime}_{i})$

##  5\. Experiments and Results

We program with python by using the third-party library PyTorch. The
experiments are conducted on a computer with an i5-8300H CPU and an Nvidia
1060 GPU. The experiment running time varies from a few minutes up to about 1
hour, depending on different agent numbers and data sizes.

###  5.1. Experiment settings

####  5.1.1. Generate Data

We randomly generate training and test data from prior distributions, and also
use the GAN network. Half of the test data is generated from the GAN and the
other half is randomly generated data. The test data size is from 10000 to
100000, and increases with $n$, as shown in Table 3. We generate new test data
for each test.

n | 4 | 5 | 6 | 7 | 8 | 9 | 10
---|---|---|---|---|---|---|---
Data size | 10000 | 20000 | 20000 | 20000 | 50000 | 50000 | 100000
Table 1. Test data size for different $n$

####  5.1.2. Batch Size

Keskar et al.(Keskar et al., 2016) show that a larger batch leads to a
dramatic degradation in the quality of the model. They investigate the cause
for this generalization drop in the large-batch regime and present numerical
evidence that supports the view that large-batch methods tend to converge to
sharp minimizers of the training and testing functions. It shows that the
large batch size converges to the sharp minimum, while the small batch size
converges to the flat.

Our experiments support Keskar’s view and we find that the loss does not
reduce with a big batch size ($\geq 1024$). We set the batch size to 64
through comparative experiments.

####  5.1.3. Order input

Arora et al. (Arora et al., 2016) shows that it is hard for ReLU to simulate
the function $max\\{a,b\\}$. A single $max\\{a,b\\}$ needs two layers and 5
nodes with exact weights and biases. So inputting sorted $\theta_{-i}$ values
to the MLP is a necessary and important step to get good results.

We order the valuations such that
$\theta_{1}\geq\theta_{2}\geq...\geq\theta_{n}$. For redistribution problems,
sorted input values will not influence strategy-proofness.

####  5.1.4. Initialization and Optimizer

We use Xavier normal initialization for the weights and $Normal(0,0.01)$ for
the bias.

We use Adam optimizer with the $learning\\_rate=0.001$ initially, and decays
by 0.98 every 100 steps by Pytorch Scheduler.

###  5.2. Results

####  5.2.1. Worst-case results

We compare our result for the worst-case optimal mechanism(GAN+MLP) with the
previously proposed mechanisms:

  * •

SBR: heuristic-based SBR mechanism (Naroditskiy et al., 2012)

  * •

ABR: heuristic-based ABR mechanism (Guo, 2016)

  * •

AMD: mechanisms derived via Automated Mechanism Design (AMD) (Guo and Shen,
2017)

  * •

AO: asymptotically optimal (AO) VCG Redistribution (Guo, 2019)

  * •

UB: the conjectured upper bounds (UB) on the efficiency ratios (Naroditskiy et
al., 2012)

n | SBR | ABR | AMD | AO | GAN + MLP | UB
---|---|---|---|---|---|---
4 | 0.354 | 0.459 | 0.600 | 0.625 | 0.634 | 0.666
5 | 0.360 | 0.402 | 0.545 | 0.600 | 0.622 | 0.714
6 | 0.394 | 0.386 | 0.497 | 0.583 | 0.592 | 0.868
7 | n too large | 0.360 | 0.465 | 0.571 | 0.626 | 0.748
8 | n too large | 0.352 | 0.444 | 0.563 | 0.654 | 0.755
9 | n too large | 0.339 | 0.422 | 0.556 | 0.682 | 0.772
10 | n too large | 0.336 | 0.405 | 0.550 | 0.623 | 0.882
Table 2. GAN+MLP Compare with previous results for worst Case

The result in Table 2 shows that our mechanism achieves better worst-case
efficiency ratios than all previous results.

The main downside of our results is that our worst-case is calculated
numerically by trying a large number of type profiles. This is a limitation
due to our neural network based approach. (Manisha et al. (Manisha et al.,
2018)’s neural network approach also evaluated the worst-case by randomly
generating a large number of type profiles. We showed that our GAN approach is
a lot more rigorous compared to simply random profile generation.) It should
be noted that SBR, ABR and AMD’s worst-cases were also calculated numerically.
(AO’s worst-case was derived analytically.) To test the stability of the
worst-case ratios, we experimented with different test sizes. Table 3 shows
that for $n=10$ with different test sizes of 10000, 20000 and 100000, $\alpha$
is stable at $0.623$.

Data size | 10000 | 20000 | 100000
---|---|---|---
$\alpha$ | 0.623 | 0.623 | 0.623
Table 3. $\alpha$ for Different test size ($n=10$)

####  5.2.2. Optimal-in-expectation results

We evaluate the expectation of $\sum_{i}{h(\theta_{-i})}/S(\theta)$ under our
mechanism obtained via MLP+FEED. As defined in our model, for a case with $n$
agents, the theoretical optimal value for $\sum_{i}{h(\theta_{-i})}/S(\theta)$
is $n-1$, so we want the expectation to be as close to $n-1$ as possible from
above (to maintain weakly budget-balance). Our results show that the obtained
mechanisms are near optimal. For example, we get $E=4.061$ vs $n-1=4$ for
$n=5$, $E=5.034$ vs $n-1=5$ for $n=6$.

Table 4 shows that for both data generated from uniform distribution and
normal distribution, the average $\sum_{i}{h(\theta_{-i})}/S(\theta)$ of our
MLP+FEED network is very close to the theoretical optimal value $n-1$, which
means that the redistribution function will return the vast majority of the
total VCG payment to the agents (particularly for a great number of agents).

MLP+FEED | Uniform(0,1) | Normal (0.5,0.1) | n-1
---|---|---|---
n=3 | 2.079 | 2.101 | 2
n=4 | 3.071 | 3.111 | 3
n=5 | 4.061 | 4.142 | 4
n=6 | 5.027 | 5.034 | 5
n=7 | 6.009 | 6.067 | 6
n=8 | 7.008 | 7.023 | 7
n=9 | 8.002 | 8.008 | 8
n=10 | 9.003 | 9.023 | 9
Table 4. $\sum_{i}{h(\theta_{-i})}/S(\theta)$ in Expectation for Different
Distributions

##  6\. Conclusion

In this paper, we consider designing optimal redistribution mechanisms for the
public project problem under two objectives: worst-case optimal and optimal-
in-expectation. With effective technical improvements on existing networks, we
train a neural network to design good redistribution functions. We use a GAN
network to generate valuation profiles to find the worst case, and feed prior
distribution into loss function to get quality gradients for the optimal-in-
expectation objective. To deal with large numbers of agents, we study
different dimension-reducing methods and supervise the network into the
existing manual mechanism as initialization. Our experiments show that for the
worst case, we could find better worst-case mechanisms compared to existing
mechanisms, and for expectation, the neural networks can derive near-optimal
redistribution mechanisms.

## References

  * (1)
  * Arora et al. (2016) Raman Arora, Amitabh Basu, Poorya Mianjy, and Anirbit Mukherjee. 2016.  Understanding deep neural networks with rectified linear units.  _arXiv preprint arXiv:1611.01491_ (2016).
  * Cavallo (2006) Ruggiero Cavallo. 2006\.  Optimal Decision-making with Minimal Waste: Strategyproof Redistribution of VCG Payments. In _Proceedings of the Fifth International Joint Conference on Autonomous Agents and Multiagent Systems_ (Hakodate, Japan) _(AAMAS ’06)_. ACM, New York, NY, USA, 882–889.  https://doi.org/10.1145/1160633.1160790
  * de Clippel et al. (2014) Geoffroy de Clippel, Victor Naroditskiy, Maria Polukarov, Amy Greenwald, and Nicholas R. Jennings. 2014.  Destroy to save.  _Games and Economic Behavior_ 86 (2014), 392–404.
  * Duetting et al. (2019) Paul Duetting, Zhe Feng, Harikrishna Narasimhan, David Parkes, and Sai Srivatsa Ravindranath. 2019\.  Optimal Auctions through Deep Learning. In _Proceedings of the 36th International Conference on Machine Learning_ _(Proceedings of Machine Learning Research, Vol. 97)_ , Kamalika Chaudhuri and Ruslan Salakhutdinov (Eds.). PMLR, Long Beach, California, USA, 1706–1715.
  * Faltings (2005) Boi Faltings. 2005\.  A Budget-Balanced, Incentive-Compatible Scheme for Social Choice. In _Lecture notes in computer science_. Springer, Berlin, 30–43.
  * Golowich et al. (2018) Noah Golowich, Harikrishna Narasimhan, and David C. Parkes. 2018\.  Deep Learning for Multi-Facility Location Mechanism Design. In _Proceedings of the Twenty-Seventh International Joint Conference on Artificial Intelligence, IJCAI-18_. International Joint Conferences on Artificial Intelligence Organization, 261–267.
  * Gujar and Narahari (2011) Sujit Gujar and Y. Narahari. 2011.  Redistribution Mechanisms for Assignment of Heterogeneous Objects.  _J. Artif. Intell. Res._ 41 (2011), 131–154.
  * Guo (2011) Mingyu Guo. 2011\.  VCG Redistribution with Gross Substitutes. In _Proceedings of the Twenty-Fifth AAAI Conference on Artificial Intelligence, AAAI 2011, San Francisco, California, USA, August 7-11, 2011_ , Wolfram Burgard and Dan Roth (Eds.). AAAI Press.  http://www.aaai.org/ocs/index.php/AAAI/AAAI11/paper/view/3733
  * Guo (2012) Mingyu Guo. 2012\.  Worst-case optimal redistribution of VCG payments in heterogeneous-item auctions with unit demand. In _International Conference on Autonomous Agents and Multiagent Systems, AAMAS 2012, Valencia, Spain, June 4-8, 2012 (3 Volumes)_ , Wiebe van der Hoek, Lin Padgham, Vincent Conitzer, and Michael Winikoff (Eds.). IFAAMAS, 745–752.  http://dl.acm.org/citation.cfm?id=2343803
  * Guo (2016) Mingyu Guo. 2016\.  Competitive VCG Redistribution Mechanism for Public Project Problem. In _PRIMA 2016: Princiles and Practice of Multi-Agent Systems - 19th International Conference, Phuket, Thailand, August 22-26, 2016, Proceedings_ _(Lecture Notes in Computer Science, Vol. 9862)_. Springer, 279–294.
  * Guo (2019) Mingyu Guo. 2019\.  An Asymptotically Optimal VCG Redistribution Mechanism for the Public Project Problem. In _Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence, IJCAI-19_. International Joint Conferences on Artificial Intelligence Organization, 315–321.
  * Guo and Conitzer (2009) Mingyu Guo and Vincent Conitzer. 2009.  Worst-case optimal redistribution of VCG payments in multi-unit auctions.  _Games and Economic Behavior_ 67, 1 (2009), 69–98.
  * Guo and Conitzer (2014) Mingyu Guo and Vincent Conitzer. 2014.  Better redistribution with inefficient allocation in multi-unit auctions.  _Artificial Intelligence_ 216 (2014), 287–308.
  * Guo and Shen (2017) Mingyu Guo and Hong Shen. 2017.  Speed up Automated Mechanism Design by Sampling Worst-Case Profiles: An Application to Competitive VCG Redistribution Mechanism for Public Project Problem. In _PRIMA 2017: Principles and Practice of Multi-Agent Systems - 20th International Conference, Nice, France, October 30 - November 3, 2017, Proceedings_ _(Lecture Notes in Computer Science, Vol. 10621)_. Springer, 127–142.
  * Keskar et al. (2016) Nitish Shirish Keskar, Dheevatsa Mudigere, Jorge Nocedal, Mikhail Smelyanskiy, and Ping Tak Peter Tang. 2016.  On large-batch training for deep learning: Generalization gap and sharp minima.  _ICLR 2017_ (2016).
  * Manisha et al. (2018) Padala Manisha, C. V. Jawahar, and Sujit Gujar. 2018\.  Learning Optimal Redistribution Mechanisms Through Neural Networks. In _Proceedings of the 17th International Conference on Autonomous Agents and MultiAgent Systems, AAMAS 2018, Stockholm, Sweden, July 10-15, 2018_ , Elisabeth André, Sven Koenig, Mehdi Dastani, and Gita Sukthankar (Eds.). International Foundation for Autonomous Agents and Multiagent Systems Richland, SC, USA / ACM, 345–353.
  * Mas-Colell et al. (1995) Andreu Mas-Colell, Michael Whinston, and Jerry R. Green. 1995\.  _Microeconomic Theory_.  Oxford University Press.
  * Moore (2006) J. Moore. 2006.  _General Equilibrium and Welfare Economics: An Introduction_.  Springer.
  * Moulin (1988) H. Moulin. 1988\.  _Axioms of Cooperative Decision Making_.  Cambridge University Press.
  * Moulin (2009) Hervé Moulin. 2009\.  Almost budget-balanced VCG mechanisms to assign multiple objects.  _JET_ 144, 1 (2009), 96–119.
  * Naroditskiy et al. (2012) Victor Naroditskiy, Mingyu Guo, Lachlan Dufton, Maria Polukarov, and Nicholas R. Jennings. 2012\.  Redistribution of VCG Payments in Public Project Problems. In _Internet and Network Economics - 8th International Workshop, WINE 2012, Liverpool, UK, December 10-12, 2012. Proceedings_ _(Lecture Notes in Computer Science, Vol. 7695)_ , Paul W. Goldberg (Ed.). Springer, 323–336.
  * Schulman et al. (2015) John Schulman, Nicolas Heess, Theophane Weber, and Pieter Abbeel. 2015.  Gradient estimation using stochastic computation graphs.  _arXiv preprint arXiv:1506.05254_ (2015).
  * Shen et al. (2019) Weiran Shen, Pingzhong Tang, and Song Zuo. 2019.  Automated Mechanism Design via Neural Networks. In _Proceedings of the 18th International Conference on Autonomous Agents and MultiAgent Systems_ (Montreal QC, Canada) _(AAMAS ’19)_. International Foundation for Autonomous Agents and Multiagent Systems, Richland, SC, 215–223.
  * Wang et al. (2021) Guanhua Wang, Runqi Guo, Yuko Sakurai, Ali Babar, and Mingyu Guo. May 2021.  Mechanism Design for Public Projects via Neural Networks. In _20th International Conference on Autonomous Agents and Multiagent Systems (AAMAS 2021, online)_.