# Traits and Transferability of Adversarial Examples against Instance
Segmentation & Object Detection

Raghav Gurbaxani Shivank Mishra
{raghavg3<EMAIL_ADDRESS>
##  1 Abstract

Despite the recent advancements in deploying neural networks for image
classification, it has been found that adversarial examples are able to fool
these models leading them to misclassify the images. Since these
classification models are now being widely deployed in self-driving cars and
drones, we provide an insight on the threat of these adversarial examples by
evaluating their characteristics and transferability to more complex models
that require Image Classification as a subtask.

We demonstrate the ineffectiveness of adversarial examples when applied to
Instance Segmentation & Object Detection. We believe that this ineffectiveness
arises from the inability of adversarial examples to withstand physical
transformations on the image such as being cropped, magnified, rotated ,or
subjected to change in lighting conditions. The loss of adversarial property
of such images while being subjected to these physical transformations raises
the question on how these adversarial examples can be put to use in the
physical world under less than perfect conditions.

Moreover, we show that there exists a small threshold below which the
adversarial property is retained while applying these input transformations.
By exceeding this threshold, the detection methods including object detectors
lose their adversarial capability.

We also extend our analysis to show that these attacks are ineffective against
Image Captioning models. Additionally, these attacks demonstrate weak cross-
network transferability across neural network architectures, _e.g._ VGG16 and
ResNet50. Lastly, we show that when an adversarial image that is generated
from one network and is passed on to another network, the resulting
adversarial image is able to fool both the networks. This may point to the
existence a universal adversarial image that is capable of fooling all
predominant network architectures.

##  2 Introduction

With recent advancements in Computer Vision, Neural Networks have been able to
achieve state-of-the-art results in Image Recognition and have been able to
outperform human-level performance. The success of Convolutional Neural
Networks [2,3,4] on the ImageNet [1] dataset (for the Large-Scale Image
Recognition Challenge) has propelled neural networks to achieve significant
results in various visual recognition tasks. Therefore, they are now widely
deployed at an unprecedented scale across a variety of applications such as
Robotics, Drones, Self-Driving Cars and Surveillance.

On the other hand, the discovery of adversarial examples [5,6] has threatened
the achievements by fooling neural networks especially during the image
classification. As pointed out by Goodfellow et. al [7] in essence, an
adversary looks to introduce a perturbation into a clean image which is almost
visually imperceptible. However, it creates a drastic change in the results of
the classification of the image. This leads to the neural network being
fooled. Despite the neural networks achieving state-of-the-art results on the
ImageNet dataset with high accuracy, the adversarial examples have been able
to fool them with a high confidence score. In order to put the grim issue in
real world perspective, an adversarial image of stop sign on road could be
misclassified as a green light.

The study of the effect of these adversarial perturbations is important,
especially if we consider that these Image detection tasks are now being
deployed widely by driverless vehicles and drones. Keeping in mind the
possible harmful ramifications of these adversaries, we introduce a
comprehensive guide on the various properties of adversarial examples and
their ability to transfer to different vision problems in order for us to
understand these attacks better.

The success of neural networks on Image Classification tasks has propelled
research in developing more complex algorithms such as Object Detection -
which classifies all the objects in the image by using Image Classifier as a
subtask and then localizes these object by forming a bounding box around them.
In recent years, researchers have also been able to perform Instance
Segmentation by further extending the Object Detection and adding a
segmentation module to it. This has allowed researchers to not only detect the
various objects in the image, but segment each pixel on those objects thus
giving us a more fine-grained classification.

With the advent of these complex models which use image classification as a
sub-task and taking into account the susceptibility of Image classifiers
against adversarial images, it becomes important to study whether the
adversarial image is able to fool these advanced pipelines as well.

In our experiments we look at the FGSM, DeepFool and Carlini & Wagner L2
attacks in a targeted and untargeted setting to study the various properties
of adversarial examples. We observed that adversarial examples lose their
capability to attack when they’re subject to cropping, magnification and
rotation. This raises the question whether there exists a threshold at which
these properties are lost. We noticed that in case of cropping and change in
lighting conditions that there exists a minimal threshold below which the
adversarial property is retained and is lost in the latter case. This
susceptibility of adversarial examples against these transformations make them
nontransferable against Object Detection & Instance Segmentation models –
which involve cropping objects in the image and classifying them individually
to their respective classes.

Our experiments on the Mask RCNN (for Instance Segmentation) and Faster RCNN &
RetinaNet (for Object Detection) corroborate our hypothesis about why they are
non-transferrable. Moreover, we also chose a third proxy task of Image
Captioning and show that adversarial examples are not transferable to them as
well. Even though Show and Fool were able to fool the captioning pipeline by
modifying the loss function, it would seem unlikely that an attacker would be
able to make such a change during an attack.

unless the attacker is able to transform the loss function of network itself
as shown by Show & Fool, which seems improbable.

In the following sections we describe our findings - We first describe the
various attack methods that we employed for testing the various models Section
3. We try out both targeted attack as well as untargeted attacking methods on
the models. Section 4 deals with the various properties of adversarial
examples and their effect when exposed to various transformations. In Section
5, we review the adversarial attacks for the task of image classification and
evaluate their transferability between different architectures. We also look
at the transferability on the complex vision models in this section. Moreover,
we also introduce a method to possibly deal with the weak transferability in
this section. In Section 6, we provide a broader outlook of the research
direction based on the reviewed literature. Finally, we draw conclusion in
Section 8.

Figure 1: Demonstrates example of adversarial misclassification

##  3 Adversarial Attacks

###  3.1 Fast Gradient Sign Method

The Fast Gradient Sign Method, as described by Goodfellow et al exploits the
linearity of neural network models in the high dimensional space. It is a form
of targeted attack – where the attacker can specify the class to which the
adversarial image should be misclassified to.

Consider _x_ as the original image, _y_ the corresponding label of _x_ ,
$\theta$ the parameters of the network and L($\theta$,_x_ ,_y_) the loss
function used to train the network.

$\eta=\epsilon sign(\nabla_{x}L(\theta,\emph{x},\emph{y}))$

The attack involves taking the derivative of the loss function in the x
direction - $\nabla_{x}L(\theta$,_x_ ,_y_) and taking the sign of the
derivative to determine the direction of the pixel change. Moreover, this is
multiplied by a small constant $\epsilon$, which ensures that the perturbation
doesn’t go too far from the original image and restricts it’s norm which gives
us the perturbation $\eta$.

This perturbation can now be added to the original image _x_ to result in the
adversarial image, capable of fooling the classifier.

$x_{adv}$ = _x_ \+ $\eta$ For our experiments we used the Cleverhans library
to generate the adversarial images on pretrained Imagenet models with with
y_target set to ’teapot’ and $\epsilon$ set to 0.3.

Figure 2: The Fast Gradient Sign Method.

(a) Clean Image

(b) Adversarial Image

Figure 3: Images: Before and after the attack

###  3.2 DeepFool

The perturbations generated by Deepfool are smaller than the perturbations
compared to FGSM in terms of norm but have similar effectiveness in terms of
fooling ratios. It also reduces the intensity of the perturbations and can
only be used in the untargeted attacking case – where the attacker can only
specify that the adversarial image is classified differently than the original
image.

This method was introduced by [5] Moosavi-Dezfooli et al. which tries to
minimize the distortion which leads to the misclassification of the
adversarial image by projecting the input image to the closest separating
hyperplane (considering each class of classifier is separated by a
hyperplane).

At every iteration, the algorithm iteratively perturbs the image by a small
vector that is calculated to take the resulting image to the boundary of the
polyhedron that is obtained by linearizing the boundaries of the region within
which the image resides. Finally, the perturbations added to the image in each
iteration are accumulated and the final perturbation changed the decision
boundary of the classifier.

For our experiments we used the Foolbox library to generate the adversarial
images on pretrained Imagenet models with with 50 iterations and $\epsilon$
set to 0.3.

Figure 4: Demonstrates example of adversarial misclassification

Although the Deepfool produces perturbations that are visibly more
indistinguishable that the FGSM, it is also a weaker attack than FGSM.

###  3.3 C&W L2

Although the distortions produced by C&W L2 are smaller than the previously
mentioned attacks, it is also a slower attack than FGSM.

Previous methods assume adversarial data can be directly fed into deep neural
networks. However, in many applications, people can only pass data through
devices (e.g., cameras, sensors). Kurakin et al. applied adversarial examples
to the physical world [20]. They extended Fast Gradient Sign method by running
a finer optimization (smaller change) for multiple iterations. In each
iteration, they clipped pixel values to avoid large change on each pixel:

where Clip x,$\xi$x,0 limits the change of the generated adversarial image in
each iteration. The adversarial examples were generated in multiple
iterations:

The authors referred to this method as Basic Iterative method. To further
attack a specific class, they chose the least-likely class of the prediction
and tried to maximize the cross-entropy loss. This method is referred to as
Iterative Least-Likely Class method:

They successfully fooled the neural network with a crafted image taken from a
cellphone camera. They also found that Fast Gradient Sign method is robust to
photo transformation, while iterative methods cannot resist photo
transformation.

Characteristics of Adversarial Attacks

It is important to study the threat of adversarial examples and how they
respond to various transformations. In this section we evaluate the
adversarial images generated from FGSM, DeepFool and C&W L2 attacks against
transformations as discussed further -

Cropping

Adversarial examples, when subjected to cropping lose their attacking
capability and switch back to their original classification as demonstrated in
the figure below. This raises two questions – 1. Whether the adversarial image
comprises of a specific global pattern and 2. Whether there exists a threshold
for cropping at which the adversarial property is lost.

For the latter question, we were able to notice that when the image was
cropped from the original size of 224 x 224 to 224 x 223 (1 pixel cropped), it
retained its adversarial property however with a significant drop in
adversarial confidence score. On further cropping, it loses its attacking
capability and the classification switches back to the original class
(Labrador retriever in this case). This also corroborates the idea that there
exists a global pattern for each adversarial image that must be preserved for
it to be able to attack classifiers.

Figure 5: Demonstrates example of adversarial misclassification

Experimental Setup
In this section we describe the datasets, the deep neural network models,
adversarial attacks methods, the evaluation metrics and the deep learning
frameworks that we used for carrying out our experiments.

a. Datasets
For generating the adversarial examples, we used the pretrained models on the
ImageNet dataset. For evaluating adversarial examples on Object detectors and
Instance Segmentation models, we used models pretrained on the MSCOCO dataset.
For evaluating Image Captioning, we used models pre- trained on the MSCOCO
dataset which comprises of 82783 training images, and 40504 and 40775
validation and test images respectively with each image corresponding to 5
captions each. We also tested our models on pretrained models on the Flikr 30k
dataset which consists of 28000 train images along with 1000 validation and
test images respectively (5 captions for each).

b. Models
For object detectors we used the famous Faster RCNN model and the RetinaNet
models, each using the ResNet 50 pipeline. We also used the Mask RCNN Instance
Segmentation module based on a ResNet 101 CNN. We used the Image Captioning
pipeline from Show and Tell based on a VGG 16 CNN. The adversarial examples
using the three attacks were generated using VGG16, ResNet 50 and Resnet 101
models pretrained on ImageNet.

c. Attacks
We employed the FGSM, DeepFool and C&W L2 using pretrained models on ImageNet.
We have mentioned our hyperparameters in Section 2. We used the FGSM and C&W
L2 attacks in a targeted setting whereas the DeepFool attack in an untargeted
setting.

d. Metrics
We used Intersection over Union (IoU) as the primary metric while evaluating
Instance Segmentation and Object Detection models. For Image Classification,
we used accuracy as a metric. For Image Captioning, we used the BLEU score.

e. Deep Learning Frameworks
We used the pretrained models on ImageNet using the Keras Library which uses
Tensorflow backend. For generating the FGSM and C&W L2 attacks, we used the
novel Cleverhans library released by Goodfellow et. al. For the DeepFool
attack we employed the Foolbox library.

Transferring attacks to Applications
Cross – Network Transferability
In this section we investigate how well do adversarial images generalize over
different neural network architectures. To do so, we use the computed
adversarial image from one model and use it to attack a different model
architecture. We tried the following experiment in a targeted attack setting.
As demonstrated in the image below we generated the adversarial attack using a
ResNet 50 based classifier which was able to fool the ResNet 50 pipeline with
high confidence. However, when we passed this adversarial image for evaluation
on a different architecture such as the VGG 16, it reverted to it’s original
(non-adversarial) class.

We achieved similar results when doing the experiments vice versa and with
different models. This indicates that the adversarial examples are weakly
transferable over neural network architectures. Object Detection

Figure 6: The Fast Gradient Sign Method.

(a) Clean Image

(b) Adversarial Image

(c) Adversarial Image

Figure 7: Images: Before and after the attack

###  3.4 Magnification

Magnification can be considered as a consequence of cropping and thus also
leads to a loss of attacking ability for the adversarial image as shown in the
figure below.

Figure 8: Maginification.

(a) Clean Image

(b) Adversarial Image

Figure 9: Instance Segmentation before and after the attack

###  3.5 Rotation

Similar to being cropped and magnified, adversarial images were not able to
preserve their attacking property on being rotated. We were interested in
finding the threshold for the susceptibility to rotation. However these images
lost their attacking capability even with a 1 degree rotation (clockwise and
anti-clockwise) as shown in the image below.

Figure 10: Rotation.

(a) Clean Image

(b) Adversarial Image

Figure 11: Instance Segmentation before and after the attack

###  3.6 Change in brightness

In case of physical adversarial examples as demonstrated by Goodfellow et. al,
they can exist in a variety of lighting conditions in the real world. Thus, we
simulate the effect of change in brightness on these images. Our experiments
divulge that adversarial examples are robust to change in exposure and retain
their attacking capability till a certain threshold. As demonstrated in the
images below- a 50% increase in exposure in the adversarial image retains the
adversarial class however a 100% increase in exposure causes it to return to
its original class.

Figure 12: The Fast Gradient Sign Method.

(a) Clean Image

(b) Adversarial Image

(c) Adversarial Image

Figure 13: Images: Before and after the attack

##  4 Experimental Setup

In this section we describe the datasets, the deep neural network models,
adversarial attacks methods, the evaluation metrics and the deep learning
frameworks that we used for carrying out our experiments.

###  4.1 Datasets

For generating the adversarial examples, we used the pretrained models on the
ImageNet dataset.

For evaluating adversarial examples on Object detectors and Instance
Segmentation models, we used models pretrained on the MSCOCO dataset. For
evaluating Image Captioning, we used models pre- trained on the MSCOCO dataset
which comprises of 82783 training images, and 40504 and 40775 validation and
test images respectively with each image corresponding to 5 captions each. We
also tested our models on pretrained models on the Flikr 30k dataset which
consists of 28000 train images along with 1000 validation and test images
respectively (5 captions for each).

###  4.2 Models

For object detectors we used the famous Faster RCNN model and the RetinaNet
models, each using the ResNet 50 pipeline. We also used the Mask RCNN Instance
Segmentation module based on a ResNet 101 CNN. We used the Image Captioning
pipeline from Show and Tell based on a VGG 16 CNN. The adversarial examples
using the three attacks were generated using VGG16, ResNet 50 and Resnet 101
models pretrained on ImageNet.

###  4.3 Attacks

We employed the FGSM, DeepFool and C&W L2 using pretrained models on ImageNet.
We have mentioned our hyperparameters in Section 2. We used the FGSM and C&W
L2 attacks in a targeted setting whereas the DeepFool attack in an untargeted
setting.

###  4.4 Metrics

We used Intersection over Union (IoU) as the primary metric while evaluating
Instance Segmentation and Object Detection models. For Image Classification,
we used accuracy as a metric. For Image Captioning, we used the BLEU score.

###  4.5 Deep Learning Frameworks

We used the pretrained models on ImageNet using the Keras Library which uses
Tensorflow backend. For generating the FGSM and C&W L2 attacks, we used the
novel Cleverhans library released by Goodfellow et. al. For the DeepFool
attack we employed the Foolbox library.

##  5 Transferring attacks to Applications

###  5.1 Cross – Network Transferability

In this section we investigate how well do adversarial images generalize over
different neural network architectures. To do so, we use the computed
adversarial image from one model and use it to attack a different model
architecture. We tried the following experiment in a targeted attack setting.

As demonstrated in the image below we generated the adversarial attack using a
ResNet 50 based classifier which was able to fool the ResNet 50 pipeline with
high confidence. However, when we passed this adversarial image for evaluation
on a different architecture such as the VGG 16, it reverted to it’s original
(non-adversarial) class.

Figure 14: Cross Network Transferability

(a) Clean Image

(b) Adversarial Image

Figure 15: Instance Segmentation before and after the attack

Resnet 50 - Labrador_retriever (58.70 $\%$)
Resnet 50 - tennis_ball (75.65$\%$)
Vgg 16 - golden_retriever (42.38$\%$)

We achieved similar results when doing the experiments vice versa and with
different models. This indicates that the adversarial examples are weakly
transferable over neural network architectures.

###  5.2 Object Detection

Figure 16: Object Detection

(a) Clean Image

(b) Adversarial Image

Figure 17: Instance Segmentation before and after the attack

###  5.3 Instance Segmentation

As an extension to the results on object detectors on the previous section, we
evaluate the attacks on the Instance Segmentation models. Despite passing the
adversarial image to the models, there is no discernable difference results of
the Mask RCNN pipelines between the original and the adversarial images. The
instance segmentation models also use the Object detection as a subtask and
adversarial attack doesn’t cause a change in result do to the susceptibility
to the cropping operation during formation of the bounding boxes.

Figure 18: Instance Segmentation.

(a) Clean Image

(b) Adversarial Image

Figure 19: Instance Segmentation before and after the attack

###  5.4 Image Captioning

In this section, we evaluate the transferability of adversarial attacks to the
Image Captioning pipeline. The goal is to recognize how the adversarial
perturbations alter the captions and their semantics. For the experimentation,
we used the captioning model with k beam search choosing k =3. Below are some
of the results obtained.

Figure 20: Demonstrates example of adversarial misclassification

Captions:
0) a man flying through the air while riding a snowboard. (p=0.021026)
1) a person jumping a snow board in the air (p=0.006679)
2) a man flying through the air while riding skis. (p=0.006655)

Figure 21: Demonstrates example of adversarial misclassification

Captions:
0) a person jumping a snow board in the air (p=0.007723)
1) a man riding a snowboard down a snow covered slope. (p=0.001448)
2) a man riding a snowboard down the side of a ramp. (p=0.000517)

Figure 22: Demonstrates example of adversarial misclassification

Captions:
0) a baseball player throwing a baseball on a field. (p=0.001463)
1) a man in a baseball uniform throwing a ball. (p=0.001192)
2) a baseball player throwing a ball on a field. (p=0.001058)

Figure 23: Demonstrates example of adversarial misclassification

Captions:
0) a baseball player swinging a bat at a ball (p=0.001510)
1) a baseball player swinging a bat on a field (p=0.000715)
2) a baseball player swinging a bat on a field. (p=0.000536)

It is evident that the introduction of adversarial perturbations does not fool
the captioning pipleline as the captions obtained by adversarial and original
images are identical. Thus, they demonstrate weak transferability.

Hypothesis: It can be inferred that the adversarial image is designed to fool
the softmax layer of the classifier. Due to the complex architecture of the
captioning model, the adversarial image may fail.

###  5.5 Attacks Capable of Fooling Multiple Networks

Extending the study described in section 5.1, we looked to explore whether
adversarial examples can fool multiple network architectures simultaneously.

To do so, we passed the adversarial image generated from one model
architecture as input for generation to the next model architecture. In some
cases, we found that the resulting adversarial image was able to fool both the
networks with a high confidence score.

In the image below, we demonstrate the aforementioned process by generating
the adversarial image using a ResNet 50 pipeline. We then passed it for
generation using a VGG16 model and found that the resulting image was able to
fool both the classification models.

However, the resulting image also suffers from more deviation from the
original and thus is visibly more altered than adversarial image meant for
only one network. This experiment would also suggest that it is possible that
an adversarial example may exist which would be capable of fooling all the
major neural network architectures.

Figure 24: Demonstrates example of adversarial misclassification

tennis_ball (99.49%) against vgg16

and tennis_ball (67.52%) against resnet 50

fools both

##  6 Future Work

###  6.1 Change Loss Function in Object Detection & Instance Segmentation

As demonstrated by Show & Fool that adversarial images can be transferred to
Image Captioning models by altering the loss function of the model. Whether
similar results can be achieved on Object Detectors and Instance Segmentation
models remains to be seen.

###  6.2 Effect of adversarial examples on Generative models

So far, adversarial images have been deployed against discriminative models.
However, it remains to be seen what the effect is of adversarial examples on
generative models such as Generative Adversarial Networks and Autoencoders.

https://arxiv.org/pdf/1801.02610.pdf

http://approximatelycorrect.com/2018/03/02/defending-adversarial-examples-
using-gans

https://arxiv.org/pdf/1710.10766.pdf

###  6.3 Effect of precision on adversarial examples

Most of the existing research on adversarial examples has been carried out on
32-bit float models. As most real-life applications use lower precision
deployments of these neural networks such as 16-bit float or 8-bit integer. It
would be intriguing to examine the effects of adversarial examples with lower
precision going down till binary.

https://openreview.net/pdf?id=HkTEFfZRb

###  6.4 Comparison of network architectures instead of fooling capability

Another possible direction would the effect of adversarial examples against
various neural network architectures. It is important to see how various
attacks fare against different architectures; which architectures are able to
defend against these attacks better.

##  7 Conclusion

We demonstrate that adversarial examples are not readily transferable to
Instance Segmentation and Object Detection models due to their sensitivity
against physical changes such as cropped, rotated or magnified. We evaluate
these tasks against the FGSM, DeepFool and C&W L2 attacks using pretrained
image classification nets on the ImageNet dataset.

In our experiments, we subjected these tasks to various input transformations
such as cropping, magnification, rotation and change in lighting conditions
and found that there exists a threshold beyond which they lose their
adversarial property (example- if an image is cropped beyond a few pixels).
When this threshold is crossed, the adversarial property is lost when
subjected to these transformations. We chose three tasks to demonstrate this
effect, namely – object detection, instance segmentation and image captioning
to show non- transferability of adversarial attacks to these tasks.

It was observed that the attacks were weakly transferable across neural
network architectures- i.e. an attack generated from VGG was unsuccessful
against fooling a ResNet based classifier (vice-versa) which corroborates the
idea that the adversarial examples are predominantly effective for the
networks from which they are generated and generally lack strong
transferability.

Lastly keeping in mind, the weak transferability, we passed an adversarial
example generated from one network (say VGG16) and passed it to generate
adversarial image on a second network (say ResNet50) and the resulting image
was able to fool both these networks. This would indicate that it would be
possible to have an adversarial image capable of fooling all the predominant
architectures.

##  8 Acknowledgement

We would like to thank the anonymous reviewer(s) for their valuable comments
and suggestions and the University of Illinois at Urbana Champaign for their
GPU computing without whose support the paper would not have been possible.

##  9 References

1\. NO Need to Worry about Adversarial Examples in Object Detection in
Autonomous Vehicles
2\. ImageNet Classification with Deep Convolutional Neural Networks.
3\. Deep Residual Learning for Image Recognition 4\. Going Deeper with
Convolutions
5\. Intriguing properties of neural networks
6\. Deep Neural Networks are Easily Fooled: High Confidence Predictions for
Unrecognizable Images
7\. Attacking Visual Language Grounding with Adversarial Examples: A Case
Study on Neural Image Captioning
8\. Cleverhans v2. 0.0: an adversarial machine learning library by N Papernot
- 2016.
9\. Explaining and Harnessing Adversarial Examples by IJ Goodfellow - 2014.
10\. DeepFool: a simple and accurate method to fool deep neural networks by SM
Moosavi-Dezfooli - 2015
11\. Universal adversarial perturbations by SM Moosavi-Dezfooli - 2016.
12\. Technical Report on the cleverhans v2.1.0 Adversarial Examples Library


# PSFs of coadded images

Rachel Mandelbaum McWilliams Center for Cosmology, Department of Physics,
Carnegie Mellon University, Pittsburgh, PA 15213, USA  Mike Jarvis Department
of Physics & Astronomy, University of Pennsylvania, 209 South 33rd Street,
Philadelphia, PA 19104-6396, USA  Robert H. Lupton Princeton University,
Princeton, NJ, USA  James Bosch Princeton University, Princeton, NJ, USA  Arun
Kannawadi Princeton University, Princeton, NJ, USA  Michael D. Murphy
McWilliams Center for Cosmology, Department of Physics, Carnegie Mellon
University, Pittsburgh, PA 15213, USA  Tianqing Zhang McWilliams Center for
Cosmology, Department of Physics, Carnegie Mellon University, Pittsburgh, PA
15213, USA  The LSST Dark Energy Science Collaboration

(August 28, 2024)

###### Abstract

We provide a detailed exploration of the connection between choice of
coaddition schemes and the point-spread function (PSF) of the resulting
coadded images. In particular, we investigate what properties of the
coaddition algorithm lead to the final coadded image having a well-defined
PSF. The key elements of this discussion are as follows:

  1. 1.

We provide an illustration of how linear coaddition schemes can produce a
coadd that lacks a well-defined PSF even for relatively simple scenarios and
choices of weight functions.

  2. 2.

We provide a more formal demonstration of the fact that a linear coadd only
has a well-defined PSF in the case that either (a) each input image has the
same PSF or (b) the coadd is produced with weights that are independent of the
signal.

  3. 3.

We discuss some reasons that two plausible nonlinear coaddition algorithms
(median and clipped-mean) fail to produce a consistent PSF profile for stars.

  4. 4.

We demonstrate that all nonlinear coaddition procedures fail to produce a
well-defined PSF for extended objects.

In the end, we conclude that, for any purpose where a well-defined PSF is
desired, one should use a linear coaddition scheme with weights that do not
correlate with the signal and are approximately uniform across typical objects
of interest.

##  1 Introduction

For many purposes, astronomers currently use and will in the future use
coadded images (i.e., images that combine the information from multiple
single-epoch observed images by registering and stacking the images in some
way) to measure the properties of stars and galaxies observed in astronomical
images. The subject of this paper is the PSF of the coadded images, and in
particular, how does the coaddition algorithm determine whether or not the
coadded image has a well-defined PSF, and how is that PSF defined? We will
henceforth use the term “coadd” to refer to the coadded image and refer to its
PSF as the “coadd PSF”. The answer to this question depends on (a) whether the
coaddition scheme is linear or not (e.g., mean versus median), and (b) what
weights are used for a linear coaddition scheme (e.g., including the object
signal in an inverse variance weight). In particular, we will show that use of
a nonlinear coaddition scheme and/or use of weights that include the signals
of individual objects when carrying out a linear coaddition scheme results in
a coadd that does not have a well-defined PSF. The authors were particularly
motivated by the challenge of analyzing data from the Vera C. Rubin
Observatory Legacy Survey of Space and Time (LSST; LSST Science Collaboration
et al., 2009; Ivezić et al., 2019), which will have hundreds of exposures at
each point within the survey footprint, but the considerations described
herein are applicable to any telescope and imaging survey.

For this purpose, we will ignore other relevant issues that occur in practice,
such as the need to obtain an accurate astrometric solution in each exposure
before producing the coadd, issues that can arise when carrying out the image
registration and resampling through interpolation111This problem is
particularly non-trivial for undersampled images, such as those from space
telescopes (see discussion in Rowe et al., 2011)., elimination of artifacts in
individual exposures, magnitude-dependent detector effects that can modify the
apparent PSF measured from stars, etc. In other words, this paper assumes that
(a) the individual images are perfectly aligned with each other, (b) the PSF
and its spatial variation are perfectly known in each individual exposure
contributing to the coadd, (c) all detector non-idealities have been perfectly
corrected, and (d) the coadd is produced at the native pixel scale. However,
we will not assume that the individual exposures have the same PSF, since the
variation in the PSFs between exposures is what makes this problem of
producing a coadd with a meaningful and well-understood PSF challenging.

This topic has been discussed in the literature in many contexts, but
unfortunately some of the key principles have not been explicitly
demonstrated. For example:

  * •

Annis et al. (2014) provide some explanation without mathematical
justification in section 3.4: “To keep track of the variance of data images,
pixel by pixel inverse variance images are a natural choice, but they produce
biases in the resulting mean. At low signal to noise the upward fluctuations
in signal are given more weight than downward fluctuations as a result of the
one-sided nature of the Poisson distribution. This bias is deterministic and
one could correct for it, but, for example, for $u$-band data with its
120$e^{-}$ of sky noise, pixels at $1\sigma$ above sky would be biased by 0.5%
and this would be a fair fraction of our photometric error budget. Another
problem is that per-pixel inverse variance weighting systematically changes
the shape of the PSF as a function of the magnitude of the object.” This is
another way of saying there is no well-defined PSF for the resulting image
when using per-pixel inverse variance weighting.

  * •

Zackay & Ofek (2017) have a thorough discussion of the impact of coaddition
choices on detection and measurements made on the coadd, but simply assert in
their introduction that the variance used for weighting should be the
“variance of all the background noise sources (e.g., background, readout
noise)”. This is an implicit rejection of weighting that includes information
about the astronomical sources such as stars and galaxies, but without
explanation. This work discusses a coaddition scheme that satisfies that
criterion and that was proposed by Nick Kaiser in an unpublished technical
note222http://pan-
starrs.ifa.hawaii.edu/project/people/kaiser/imageprocessing/im++.pdf.

  * •

Mandelbaum (2018) asserts without explanation in section 2.4 that depending on
how coaddition is carried out, “the coadd may not even have a well-defined PSF
at each point (for example, inverse variance weighting that depends on the
total flux, or use of a median for the coadd). The answer to this challenge is
to generate the coadd in a principled way that results in a single well-
defined PSF at each point.” This statement is consistent with the above two
references, but is provided without explanation.

  * •

Bosch et al. (2018) asserts without illustration in section 3.3.2 that “For
PSF model coaddition to be valid, the operation used to combine all input
pixels at each point on the coadd image must be strictly linear – robust
estimators such as the median or sigma-clipped means cannot be used. Nonlinear
estimators do not just prevent PSF coaddition from working, however; they
prevent the coadded image from even having a well-defined coadd PSF. Any
estimator that rejects individual pixel outliers will tend to reject pixels in
the cores of stars on the best-seeing exposures, and brighter stars will
experience more rejection, giving them a different profile than fainter stars.
It should be emphasized that this occurs even in the absence of noise, and
even with extremely weak outlier rejection (e.g. clipping at $10\sigma$). All
robust estimators start from the ansatz that all input values are drawn from
the same underlying distribution, and convolution with different PSFs means
that they are not.” This is an even more general statement, rejecting all
nonlinear coaddition schemes; an explanation is given only for the problems
with outlier rejection during coaddition, and not for the more general
statements at the start of the quoted material.

  * •

A technical note about coaddition algorithms333https://dmtn-015.lsst.io/ by
Jim Bosch (for Rubin Observatory Data Management) states that when producing a
coadd via a weighted summation over individual exposures, “It is equally
important that the weight function vary only on scales much larger than the
size of the PSF; without this the coadd has no well-defined effective PSF.”

  * •

Rowe et al. (2011) note in section 2.2 that they derive a linear coadd
estimator because this “ensures that the final PSF … is independent of image
brightness, so the stars that were not used in the PSF determination are
tracers of how well the image synthesis worked.”

In this paper we aim to provide some of the mathematical justification and
physical intuition behind these statements, beginning with defining formalism
in Section 2, then discussing linear and nonlinear coaddition schemes in
Sections 3 and 4, respectively.

##  2 Formalism

The most basic definition of a “point-spread function” is that it is a
function that describes the way that a bundle of rays, which would be incident
on a single point, say $(x_{0},y_{0})$, spreads out spatially around that
point. That is, if the incident profile is a Dirac delta function
$\delta(x-x_{0},y-y_{0})$, then the function $P(x,y)$ describes the resulting
distribution of rays relative to the nominal point of incidence:

$\delta(x-x_{0},y-y_{0})\rightarrow P(x-x_{0},y-y_{0}).$ (1)

Here the $(x,y)$ coordinates are not necessarily detector coordinates; they
could instead be some idealized sky coordinate system after correction for
astrometric distortions. Across a typical astronomical image, this function
$P$ may vary as a function of position $x_{0},y_{0}$, but it typically does so
fairly slowly. Under the approximation that the functional form of the PSF is
essentially unchanging across the extent of a single galaxy or other object or
area of interest, the observed image can be described as a convolution:

$\displaystyle I(x,y)$ $\displaystyle=T(x,y)\otimes P(x,y),$ (2)
$\displaystyle\equiv\int_{-\infty}^{\infty}\mathrm{d}x^{\prime}\int_{-\infty}^{\infty}\mathrm{d}y^{\prime}T(x^{\prime},y^{\prime})P(x-x^{\prime},y-y^{\prime})$
(3)

where $\otimes$ denotes convolution and $T(x,y)$ is the underlying true
surface brightness distribution of the relevant sky scene.

###  2.1 The Coadd PSF

In this work we assume that multiple background-subtracted images $I_{i}(x,y)$
(with subscript $i$ indexing the individual images) are being combined to
produce a coadded image $I_{\text{coadd}}(x,y)$. We also presume that each
image has the same true underlying sky scene444This means that technically,
our results only apply to static sources. Following the results of this paper,
it seems that the coadd PSF for time-varying sources will differ depending on
the light curve (e.g., in the limit that a variable or transient source is
very bright in just one exposure and faint in all the others, its coadd PSF
will essentially be that of the exposure where it is very bright). This effect
could even cause more basic issues for LSST. The initial survey plan was that
each visit (observation at a particular point) would consist of two 15-second
“snaps” which would be combined into a single 30-second exposure. It is still
to be decided whether this will be done, but if it is, then the issue of the
PSF for time-varying objects in a coadd is even relevant to how the two snaps
are combined. Exploring these issues for time-varying sources is beyond the
scope of this paper. $T(x,y)$ (hence it has no index $i$) but a potentially
different effective PSF (c.f. Sec 2.2) $P_{i}(x,y)$. Finally, as mentioned
above, we assume that the form of the PSF may vary on large spatial scales,
but is assumed to not change significantly across the extent of an individual
star or galaxy or other local area of interest in the image. Within a portion
of the image small enough that these assumptions hold, the image amplitude
prior to the addition of noise may be represented as

$I_{i}(x,y)=T(x,y)\otimes P_{i}(x,y).$ (4)

We investigate the following question: under what circumstances does the
coaddition process of combining multiple images $I_{i}(x,y)$ produce a coadded
image $I_{\text{coadd}}(x,y)$ such that there is a well defined function
$P_{\text{coadd}}(x,y)$ that satisfies the condition

$I_{\text{coadd}}(x,y)=T(x,y)\otimes P_{\text{coadd}}(x,y),$ (5)

or its Fourier-space equivalent

$\widetilde{I}_{\text{coadd}}(k_{x},k_{y})=\widetilde{T}(k_{x},k_{y})\widetilde{P}_{\text{coadd}}(k_{x},k_{y})?$
(6)

As we will see, the existence of a function $P_{\text{coadd}}(x,y)$ satisfying
this condition puts quite strong requirements on the nature of the coaddition.
For instance, it implies that any stellar (point-source) images must scale
linearly with the flux of the delta-function source. This requirement is the
origin of several tests we will carry out in Sections 3 and 4. It also imposes
requirements on extended sources, which we will explore further in Section
4.3.

Each image is presumed to have a spatially-varying background level
$B_{i}(x,y)$ that has been subtracted prior to PSF measurement and coaddition.
This background level may potentially be used to define a spatially-varying
weight function for coaddition. Moreover, while not included explicitly as a
separate term, each image is presumed to have noise (the specific form of
which is irrelevant to the rest of this discussion).

###  2.2 Sampling

The above formalism treats the image as a continuous function of $(x,y)$. In
reality, one deals with pixellated images $I_{i}[m,n]$ where
$m,n\in\mathbb{Z}$ denote the pixel index. A pixellated image is not merely
the incident flux density sampled at discrete $(x,y)$ locations, i.e.,
$I_{i}[m,n]\neq I_{i}(x=m,y=n)$. It is, rather, an integral of this function
within the pixel area. For simplicity, we illustrate this with a 1-dimensional
flux density $I_{i}(x)$, but the argument is applicable to 2-dimensional
images as well. For a uniform pixel spacing $a$, the pixel value at the
$n^{\text{th}}$ pixel $I_{i}[n]$ is given by

$I_{i}[n]=\int_{(n-0.5)a}^{(n+0.5)a}\mathrm{d}x^{\prime\prime}I_{i}(x^{\prime\prime}).$
(7)

As a consequence, the total flux is conserved since

$\sum_{n}I_{i}[n]=\int_{-\infty}^{\infty}dx^{\prime\prime}I_{i}(x^{\prime\prime}).$
(8)

We now show that the pixellated image is still a convolution of the underlying
scene with an effective PSF, sampled at the centers of pixels (see also Lauer,
1999; Rowe et al., 2015). Starting from Eq. (7),

$\displaystyle I_{i}[n]$
$\displaystyle=\int_{(n-0.5)a}^{(n+0.5)a}\mathrm{d}x^{\prime\prime}I_{i}(x^{\prime\prime})$
$\displaystyle=\int_{(n-0.5)a}^{(n+0.5)a}\mathrm{d}x^{\prime\prime}\int_{-\infty}^{\infty}\mathrm{d}x^{\prime}T(x^{\prime})P_{i}(x^{\prime\prime}-x^{\prime}).$

As before, $T$ is the underlying light profile of the sky scene, and $P_{i}$
is the PSF in image $i$. Next we explicitly introduce the pixel response
function, $\Pi(u)=1$ if $|u|<a/2$ and 0 otherwise:

$I_{i}[n]=\int_{-\infty}^{\infty}\mathrm{d}x^{\prime}T(x^{\prime})\left[\int_{-\infty}^{\infty}\mathrm{d}x^{\prime\prime}P_{i}(x^{\prime\prime}-x^{\prime})\Pi(x-x^{\prime\prime})\right|_{x=na}.$

Finally, we change the variable of integration from $x^{\prime\prime}$ to
$q=x^{\prime\prime}-x^{\prime}$ in the inner integral,

$\displaystyle I_{i}[n]$
$\displaystyle=\int_{-\infty}^{\infty}\mathrm{d}x^{\prime}T(x^{\prime})\left[\int_{-\infty}^{\infty}\mathrm{d}qP_{i}(q)\Pi\left((x-x^{\prime})-q\right)\right|_{x=na},$
$\displaystyle=\int_{-\infty}^{\infty}\mathrm{d}x^{\prime}T(x^{\prime})P_{\text{eff},i}(x-x^{\prime})\big{|}_{x=na},$
(9)

where we define the effective PSF as the convolution of the PSF with the pixel
response:

$P_{\text{eff},i}(x)\equiv P_{i}(x)\otimes\Pi(x).$ (10)

We can thus consider the discrete pixel values as samples of the underlying
scene convolved with the effective PSF, evaluated at the pixel centers.
However, for simplicity of notation, in the rest of this note we simply write
$P(x,y)$ rather than $P_{\text{eff}}(x,y)$ to denote the effective PSF.

Since the convolution operation is associative, it is possible to include
sensor effects such as charge diffusion, interpixel capacitance etc. as part
of the effective PSF, as long as those effects act as convolutions (which is
not the case for all sensor effects). We use the phrase “effective PSF” to
refer to the convolution of the PSF from the atmosphere and telescope optics
with the pixel response function and with sensor effects that act as
convolutions. We prefer this convention over the alternative, i.e. not
including the pixel response and then explicitly integrating over the pixel
response function, as it makes the bookkeeping cleaner for coadded or
resampled images.

Once we redefine $P(x,y)$ to really mean the effective PSF
$P_{\text{eff}}(x,y)$, then $I_{i}(x,y)$ as defined in Eq. (5) is a continuous
function that, when sampled at the pixel centers, gives the observed image
values. Throughout the rest of the paper, we neglect this final step of
sampling this function at the pixel centers, as it is irrelevant to the points
we are making. The results that we find apply to $I(x,y)$ will apply equally
well to the real $I[m,n]$ data, regardless of whether the PSF is Nyquist
sampled or other similar details.

##  3 Linear coaddition schemes

In this section, we presume that the images are summed with weight functions
$w_{i}(x,y)$, such that

$I_{\text{coadd}}(x,y)=\frac{\sum_{i}w_{i}(x,y)I_{i}(x,y)}{\sum_{i}w_{i}(x,y)}.$
(11)

The denominator ensures flux conservation in the coaddition process by
enforcing that the relative weights $w_{i}/\sum_{i}w_{i}$ always sum to 1.
Note that this is not the most general linear coadd estimator; a matrix
product would be more general, and in the case of undersampled images, a
matrix product555Here the “matrix product” would connect pixels in the coadd
with pixels in the input image. For a list of $n_{\text{in}}$ input pixel
values $\mathbf{d}_{\text{in}}$ in the set of images to be coadded, and a list
of $n_{\text{out}}$ desired coadd pixel values $\mathbf{d}_{\text{out}}$, the
most general relationship is
$\mathbf{d}_{\text{out}}=\mathbf{W}\mathbf{d}_{\text{in}}$, where $\mathbf{W}$
is an $n_{\text{out}}\times n_{\text{in}}$ matrix. In practice, even for
approaches employing a matrix product, a variety of approximations are made to
this general (but computationally unwieldy) formulation. is actually necessary
to get a well-defined coadd PSF (Rowe et al., 2011). However, for simplicity
of notation we use Eq. (11) throughout this paper.

We begin with a direct illustration of the impact of weighting choices in a
linear coaddition scheme, and then follow with a more formal demonstration of
the circumstances under which the coadd has a well-defined PSF.

###  3.1 Direct illustration

Figure 1: Result of coadding three one-dimensional images of stars with
different PSFs using per-pixel inverse variance weighting that depends on the
signal, for two different choices of the star flux $f_{\text{star}}$. The
figure shows the three input 1D images in blue, orange, and green for the
three star images $I_{1}(x)$, $I_{2}(x)$, and $I_{3}(x)$ (Gaussians with
$\sigma=0.5$, $1.5$ and $3.0$, respectively), all normalized to integrate to
one for convenience of comparison. The resulting coadded images
$I_{\text{coadd}}(x)$ for two different stellar fluxes are also shown, each
produced assuming a constant background level of $10$ in each image, with
inverse variance weighting from Eq. (12) including both background and signal.
As shown, the normalized coadded images differ depending on whether the images
are background-dominated (lower star flux, shown in black) or signal-dominated
at the star center (higher star flux, shown in pink).

A natural choice for weighting is inverse variance weighting, to achieve an
optimal signal-to-noise in the detection of objects in the coadd666Inverse
variance weighting is statistically optimal unless one adopts more general
linear coadd estimators involving a matrix product, as in the Kaiser
approach.. We therefore define a fully spatially varying optimal weight for
each image:

$w_{i,\text{opt}}(x,y)=\frac{1}{B_{i}(x,y)+I_{i}(x,y)},$ (12)

where $B_{i}$ is the background level as before. Combining Eqs. (11) and (12),
the coadded image with this optimal weight is defined as

$I_{\text{coadd,opt}}(x,y)=\frac{\sum_{i}I_{i}(x,y)/\left[B_{i}(x,y)+I_{i}(x,y)\right]}{\sum_{i}1/\left[B_{i}(x,y)+I_{i}(x,y)\right]}.$
(13)

To illustrate the impact of weight functions, we presume we are viewing an
image of a single star. In this case, we can assume that the sky scene
consists of a Dirac delta function at the center of our image,
$(x_{0},y_{0})=(0,0)$, so $T(x,y)=f_{\text{star}}\delta(x,y)$, and therefore
$I_{i}(x,y)=f_{\text{star}}P_{i}(x,y)$ where $f_{\text{star}}$ is the total
star flux. We further assume that the PSFs and backgrounds in the images may
in general differ from each other, but are spatially constant, i.e.,
$B_{i}(x,y)=B_{i}$. For simplicity, we consider three PSFs that are described
as Gaussians with $\sigma=0.5$, $1.5$ and $3.0$ (arbitrary units), and using
the same background level of $B=10$ for each. We then ask the question of what
does the coadd look like for different values of $f_{\text{star}}$? The result
is shown in Figure 1.

As shown, the resulting coadded image differs depending on whether the star
flux is such that the images are background-dominated at all positions, or
whether there are some positions near the center of the star that are signal-
dominated. This is a key result: in reality, there are objects of different
fluxes throughout the images, and this figure suggests that once the images
have different PSFs, the dependence of the weights on the sources’ fluxes
means that objects will appear differently in the coadds depending on their
brightness, even in the very simple case that we only have three images with
the same constant background level and with a PSF that does not vary spatially
within the images. Once that is the case, there is no well-defined PSF for
these images, as defined in Eq. (5). Therefore, despite the statistical
optimality of this weighting scheme, it is fundamentally flawed from the
standpoint of systematics; the lack of a well-defined coadd PSF will induce
systematic biases in the characterization of object properties.

Ultimately the origin of these issues is that the weights that define the
coadd in Eq. (13) depend on the input signal in the objects in the sky scene.
It seems possible that weights that depend only on quantities that are
unrelated to the sky scene (example: background level and average seeing size)
will not induce the same problems, which is the origin of the choices made in
Annis et al. (2014) and Zackay & Ofek (2017). This statement that weighting by
the background level should not cause problems with the coadd PSF is true
under the assumption that the background level is dominated by other factors
besides the astronomical objects (stars and galaxies) of interest, e.g.,
scattered light in the Earth’s atmosphere or zodiacal light, which is
generally the case. By extension, a coadd with weights that depend on object
signals will come closer to having a well-defined PSF in the limit that all
objects are background-dominated.

###  3.2 Formal demonstration

A more principled mathematical view of the issue comes from the relationship
between these quantities in Fourier space. Starting from Eq. (11), we can use

$\displaystyle I_{\text{coadd}}(x,y)$
$\displaystyle=\frac{\sum_{i}w_{i}(x,y)I_{i}(x,y)}{\sum_{i}w_{i}(x,y)}$
$\displaystyle=\frac{\sum_{i}w_{i}(x,y)\text{InvFT}[\widetilde{I}_{i}(k_{x},k_{y})]}{\sum_{i}w_{i}(x,y)}$
$\displaystyle=\frac{\sum_{i}w_{i}(x,y)\text{InvFT}[\widetilde{T}(k_{x},k_{y})\widetilde{P}_{i}(k_{x},k_{y})]}{\sum_{i}w_{i}(x,y)}.$

And Eq. (5) tells us that we want to know when is there some function
$\widetilde{P}_{\text{coadd}}(k_{x},k_{y})$ such that

$I_{\text{coadd}}(x,y)=\text{InvFT}[\widetilde{T}(k_{x},k_{y})\widetilde{P}_{\text{coadd}}(k_{x},k_{y})]?$

We could imagine trying to find $\widetilde{P}_{\text{coadd}}$ via division,

$\widetilde{P}_{\text{coadd}}(k_{x},k_{y})=\frac{\widetilde{I}_{\text{coadd}}(k_{x},k_{y})}{\widetilde{T}(k_{x},k_{y})}.$
(14)

However, for there to be a valid PSF, it needs to be the case that the
(Fourier-transformed) true scene $\widetilde{T}$ divides out of this ratio. In
other words, the ratio should be independent of $\widetilde{T}$. It seems that
for a completely general set of weights to result in a well-defined PSF for
the coadd, $\widetilde{I}_{\text{coadd}}(k_{x},k_{y})$ must include the scene
in the form of a straightforward product. Or, equivalently, it must be
possible to write $I_{\text{coadd}}(x,y)$ as the convolution of the true scene
and some other function. For a linear coaddition scheme, this will only occur
under the following conditions:

  * •

Same PSF in each image: If all the component exposures have the same PSF
(which could be approximately true, e.g., for space-based observations, where
there is no atmosphere, if there is very little dithering or other effects
such as thermal “breathing” between the exposures), then we have
$P_{i}(x,y)=P(x,y)$. In this case, Eq. (11) reduces to

$\displaystyle I_{\text{coadd}}(x,y)$
$\displaystyle=\frac{\sum_{i}w_{i}(x,y)T(x,y)\otimes
P(x,y)}{\sum_{i}w_{i}(x,y)}$ $\displaystyle=T(x,y)\otimes P(x,y).$

Since each image has the same PSF, they are all the same (modulo noise). Hence
the sum over the spatially varying weight function in each image cancels out
of the numerator and denominator, and has no impact on the PSF of the
resulting coadd. The coadded image has a well-defined PSF: it is the same as
the PSF in each input image.

  * •

Spatially constant weight functions: If the weight for each exposure is not a
function of position on the exposure, then $w_{i}(x,y)=w_{i}$. In this case,
we can rewrite Eq. (11) as

$\displaystyle I_{\text{coadd}}(x,y)$
$\displaystyle=\frac{\sum_{i}w_{i}I_{i}(x,y)}{\sum_{i}w_{i}}$
$\displaystyle=\frac{\sum_{i}w_{i}\left[T(x,y)\otimes
P_{i}(x,y)\right]}{\sum_{i}w_{i}}$
$\displaystyle=T(x,y)\otimes\frac{\sum_{i}w_{i}P_{i}(x,y)}{\sum_{i}w_{i}},$
(15)

where the last step relies on the linearity of convolution. Eq. (15) shows
that in this case, the scene is included in $I_{\text{coadd}}(x,y)$ as a
direct convolution, which means the coadd has a well-defined PSF,

$P_{\text{coadd}}(x,y)=\frac{\sum_{i}w_{i}P_{i}(x,y)}{\sum_{i}w_{i}}.$ (16)

An example of a valid weighting scheme of this type would be if the weights
were set to a spatially constant value determined based on the average PSF
size in each contributing exposure.

  * •

Weights that are independent of the true scene: If the weights are based
solely on imaging characteristics (e.g., the background level, electronic read
noise, etc.), and not anything correlated with the true scene, then the
weights $w_{i}(x,y)$ are uncorrelated with $T(x,y)$. In this case, we have
$\langle(w_{i}(x,y)-\bar{w}_{i})T(x,y)\rangle=0$, where the expectation value
$\langle...\rangle$ averages over the spatial extent of a large survey and
$\bar{w}_{i}=\langle w_{i}(x,y)\rangle$ is the spatial average of the weights.

Define $P_{\text{coadd}}$ to be the weighted average PSF:

$P_{\text{coadd}}(x,y)=\frac{\sum_{i}w_{i}(x,y)P_{i}(x,y)}{\sum_{i}w_{i}(x,y)}.$
(17)

It will not be the case now that we can satisfy Eq. (5) specifically at every
location, since $\widetilde{T}$ will not cancel out in Eq. (14), but we can
satisfy it in an expectation-value sense, averaging over all locations in a
survey (calculation shown in 1D rather than 2D for notational simplicity):

$\displaystyle\left\langle I_{\text{coadd}}(x)-T(x)\otimes
P_{\text{coadd}}(x)\right\rangle$
$\displaystyle=\left\langle\frac{\sum_{i}w_{i}(x)\int
dx^{\prime}T(x^{\prime})P_{i}(x-x^{\prime})}{\sum_{i}w_{i}(x)}-\int\\!dx^{\prime}\frac{T(x^{\prime})\sum_{i}w_{i}(x-x^{\prime})P_{i}(x-x^{\prime})}{\sum_{i}w_{i}(x-x^{\prime})}\right\rangle$
$\displaystyle=\left\langle\sum_{i}\int\\!dx^{\prime}T(x^{\prime})P_{i}(x-x^{\prime})\left(\frac{w_{i}(x)}{\sum_{i}w_{i}(x)}-\frac{w_{i}(x-x^{\prime})}{\sum_{i}w_{i}(x-x^{\prime})}\right)\right\rangle$
$\displaystyle=\sum_{i}\int\\!dx^{\prime}T(x^{\prime})P_{i}(x-x^{\prime})\left(\left\langle\frac{w_{i}(x)}{\sum_{i}w_{i}(x)}\right\rangle-\left\langle\frac{w_{i}(x-x^{\prime})}{\sum_{i}w_{i}(x-x^{\prime})}\right\rangle\right)$
$\displaystyle=0$

where the second step (passing $T$ and $P_{i}$ out of the expectation
brackets) relies on $w_{i}$ being uncorrelated with $T$ and $P_{i}$777Note
that technically in addition to being uncorrelated with $T(x,y)$, the weights
must also be uncorrelated with the PSF profile values $P_{i}(x,y)$, so
$\langle(w_{i}(x,y)-\bar{w}_{i})P_{i}(x,y)\rangle=0$. For example, weight
functions that depend on the spatial profile of the PSF around some particular
central point would be disallowed. As this would be a bizarre choice for
coaddition weights, we ignore this fairly pedantic technical point for the
remainder of the paper. We note in particular that weights that depend on
overall PSF characteristics such as size or shape are acceptable..

Thus, in this case, we have the slightly modified condition:

$\left\langle I_{\text{coadd}}(x,y)\right\rangle=\left\langle T(x,y)\otimes
P_{\text{coadd}}(x,y)\right\rangle$ (18)

Each location on the image will not specifically follow Eq. (5), but using
this equation as a model of the observed coadd image will be unbiased when
applied to a large ensemble of sources. For many purposes, this is a
sufficient definition for there being a well-defined PSF.

These are the only options we were able to identify that guarantee that the
linear coadd defined in Eq. (11) can be rewritten as the convolution of the
true sky scene with some other function. We do not believe that there are
others where it is possible to guarantee that the coadded image will have a
well-defined PSF888Note that one could imagine devising spatially varying
weights that depend on $T(x)$, yet cause the various terms in $\langle
I_{\text{coadd}}(x)-T(x)\otimes P_{\text{coadd}}(x)\rangle$ to cancel for a
specific choice of $T(x)$. However, they would not be guaranteed to cancel in
the real-world case that $T(x)$ is not perfectly known, unless the above-
mentioned conditions apply..

For real data, it will never be the case that each image has precisely the
same PSF; therefore, selecting spatially constant weight functions, or at
least weight functions that are independent of the signal $T(x,y)$, provides a
viable route to ensuring the coadd has a well-defined PSF. In practice, the
last option seems to be the most common, since the various references in
Section 1 of this paper accept the idea of weighting by the spatially-varying
background, with the Rubin Data Management technical note providing an
approximate criterion: the weight function should “vary only on scales much
larger than the size of the PSF.” This criterion helps guarantee that the
weights will be uncorrelated with the sources, at least for the majority of
the astronomical objects in the images, with spatial extents that are
typically just a few times the PSF size (for example, the small galaxies that
are typically used for dark energy science).

##  4 Nonlinear coaddition schemes

There are several possible nonlinear coaddition schemes that one might
consider adopting in reality. While the discussion in Sec. 3.2 makes it
challenging to envision a nonlinear approach that would result in a well-
defined PSF for the coadd, we discuss the specific issues with two potential
nonlinear approaches below before explaining mathematically why it is
impossible for a nonlinear coadd scheme to have a well-defined PSF.

###  4.1 Median coadds

In order to avoid contamination by image artifacts, one possible “robust
statistic” is a median rather than a mean. In Figure 2 we illustrate the
hazards of a median approach. We show the images for three one-dimensional
Gaussian images with different PSFs, along with the profile generated by a
median coadd. We generated one million noisy realizations of the images with a
sky background of 100, used them to make coadds, and then took the average of
that coadd across all of those realizations to beat down the noise. This was
done for both low- and high-flux stars (flux = $1$ and $10^{5}$ respectively).
As shown, the high-flux case results in a star image with discontinuous slope,
as expected from the medians of the three noise-free star images. In contrast,
in the low-flux case, the star image resembles an unweighted mean coadd. This
illustrates that the PSF depends on the object magnitude, meaning that the
coadd does not have a single well-defined PSF.

Figure 2: Result of applying a median coaddition procedure to three one-
dimensional images of stars with different PSFs. First, we show the images for
three one-dimensional Gaussian images with different PSFs (denoted $I_{i}(x)$
in the legend) – with the flux divided out to avoid stars of different fluxes
requiring very different axes. The other curves show average coadd profiles
from many noisy realizations, including a sky background of 100. The black
dashed curve shows an unweighted mean coadd profile (for any star flux), which
has a well-defined PSF. The solid grey curve shows the result for a median
coadd for a low-flux (background-dominated) star, which as shown looks very
similar to the mean coadd. Finally, the pink curve shows the result for the
median coadd for a high-flux (signal-dominated) star. The median coadd curves
are shown as very broad but semi-transparent lines in order to make it more
clear that for various $x$ values they overlap various of the single-exposure
images $I_{i}(x)$ or (for the low-flux case) the $\langle
I_{\text{mean}}(x)\rangle$ curve. The fact that the pink and grey curves
differ demonstrates that there is no well-defined PSF for a median coadd.
Figure 3:  Visualization of biases in PSF photometry for the median coaddition
scheme, using 3 differing input exposures for the coadd. Different colors
indicate a different PSF size distribution. The legend shows a factor
indicating the ratio of worst to best seeing size, along with the 3 PSF FWHM
values (Full-Width at Half-Maximum in arcsec, where the pixel scale is 0.2
arcsec). An image of a star that is very bright compared to the background
level is used as the PSF model, $P(x)$, and fit to the image of a faint star
to measure the PSF flux, as described in text. The curves show the average
fractional bias in the estimated faint star flux, using 1000 independent noise
realizations, as a function of the ratio of faint-to-bright star flux ratio.
The plot is shown on a symlog y-axis scale; light dotted grey lines indicate
the range within which the y axis is linear, while the dark dotted grey line
shows the ideal (bias-free) result of 0. As expected, all curves converge to 0
as the faint star flux increases to match that of the bright star, but
significant biases are present at lower flux, especially when the distribution
of PSF sizes contributing to the coadd is relatively broad.

For a quantitative analysis of the impact of this coaddition scheme with more
realistic inputs, Figure 3 shows the bias introduced by the use of a median
coadd when measuring point source photometry using PSF model flux
measurements. This algorithm uses a normalized bright star image as a PSF
model, and finds the flux that best describes the light profiles of other
stars, minimizing the $\chi^{2}$ as a summation over image pixels with the
same weight assigned to each pixel:

$\chi^{2}=\sum_{\text{pixels}}\frac{(I_{\text{median}}(x)-f_{*}P(x))^{2}}{\sigma^{2}},$
(19)

where $I_{\text{median}}(x)$ is the image of the star for which we wish to
measure the flux; $f_{*}$ is the fitting parameter, representing the apparent
flux; $P(x)$ is the PSF model evaluated using a very bright star image; and
$\sigma$ is the background noise level in the image, which we take to be a
constant. Figure 3 shows the fractional bias in these PSF flux estimates, as a
function of the faint-to-bright star flux ratio, with means and uncertainties
shown using fits to 1000 independent noise realizations. As expected, the bias
is zero when the “faint” star flux approaches the bright star flux; but for
fainter star fluxes, biases in photometry appear, with a magnitude that
depends on the distribution of seeing sizes in the images contributing to the
coadd. When the seeing varies by only 20% among the different images, biases
in the PSF photometry are $\lesssim$1%, comparable to the uncertainties. As
the range of PSF sizes increases, so does the bias in PSF photometry,
eventually reaching $\sim$10% when the range is a factor of 2–3. Like other
ground-based surveys, LSST is likely to have ratio around 2, for which the
median coaddition scheme would cause unacceptable levels of bias in point-
source photometry.

###  4.2 Clipped-mean coadd

Another obvious “robust statistic” to try is a clipped mean, i.e., defining
some $\sigma$ value (through an inter-quartile range or direct calculation),
and omitting pixels from an image if they deviate from the mean by more than
$N\sigma$, possibly working iteratively until convergence is achieved. This is
often a robust method for the problem of removing outliers from data with
occasional glitches, e.g. cosmic rays in astronomical data.

However, as noted in Sec. 1, Bosch et al. (2018) has pointed out the issues
with clipped means. Functionally, they are the solution to a different
problem, one where all of the data are drawn from the same underlying
distribution. In the case of coaddition, different images have different PSFs,
so each observed image $I_{i}(x,y)$ is drawn from a different distribution;
therefore, a clipped mean is not actually a robust statistic for the problem
of removing outlier data in the input images. For example, images with a
smaller PSF will preferentially have higher pixel values at the centers of
bright stars and galaxies, which will preferentially fail such an $N\sigma$
criterion orders of magnitude more frequently than one would expect for a
Gaussian distribution in cases where the data are drawn from the same
underlying distribution. Omitting data preferentially as a function of PSF
size, especially doing so differently for different parts of an object (i.e.,
core vs. wings), will naturally result in a coadd with no well-defined PSF.
This is true even for relatively large values such as $N=10$.

###  4.3 Extended Sources

The above arguments about median and clipped-mean coadds seem to imply that
they are only a problem at moderately low signal-to-noise. At high signal-to-
noise, the profiles of stars are nearly independent of the flux of the star
for both median and clipped-mean coadds. However, there is an additional
problem with nonlinear coaddition schemes when considering extended sources
(i.e., galaxies), even when limiting to high signal-to-noise objects.

We typically determine the PSF profile by looking at the images of stars,
since these represent the response of the imaging process to a delta function.
This essentially acts like a Green’s function of the imaging “operator”.

$P_{i}(x,y)=\mathbf{F}_{i}\delta(x,y),$ (20)

where $\mathbf{F}_{i}$ is an operator representing some potentially
complicated imaging process.

We can then model the same process acting on an extended object

$I_{i}(x,y)=\mathbf{F}_{i}T(x,y)$ (21)

as a convolution

$\displaystyle I_{i}(x,y)$ $\displaystyle=\mathbf{F}_{i}\int
dx^{\prime}dy^{\prime}T(x^{\prime},y^{\prime})\delta(x-x^{\prime},y-y^{\prime})$
$\displaystyle=\int
dx^{\prime}dy^{\prime}T(x^{\prime},y^{\prime})\mathbf{F}_{i}\delta(x-x^{\prime},y-y^{\prime})$
$\displaystyle=\int
dx^{\prime}dy^{\prime}T(x^{\prime},y^{\prime})P_{i}(x-x^{\prime},y-y^{\prime})$
$\displaystyle=T(x,y)\otimes P_{i}(x,y),$ (22)

which is valid so long as $\mathbf{F}_{i}$ is a linear operator.

When we coadd multiple images $I_{i}(x,y)$ into $I_{\mathrm{coadd}}(x,y)$,
this coadding process essentially becomes part of the imaging operator
$\mathbf{F}_{\mathrm{coadd}}$:

$\displaystyle I_{\mathrm{coadd}}(x,y)$
$\displaystyle=\mathrm{Coadd}\left(\\{I_{i}(x,y)\\}\right)$
$\displaystyle=\mathrm{Coadd}\left(\\{\mathbf{F}_{i}(x,y)T(x,y)\\}\right)$
$\displaystyle\equiv\mathbf{F}_{\mathrm{coadd}}T(x,y).$ (23)

When applied to stars, after normalizing the images by their fluxes, we can
call the resulting profile $P_{\mathrm{coadd}}$, but as we will see, this is
not really a PSF when $\mathbf{F}_{\mathrm{coadd}}$ is nonlinear.

$\displaystyle P_{\mathrm{coadd}}(x,y)$
$\displaystyle=\mathrm{Coadd}\left(\\{P_{i}(x,y)\\}\right)$
$\displaystyle=\mathrm{Coadd}\left(\\{\mathbf{F}_{i}(x,y)\delta(x,y)\\}\right)$
$\displaystyle=\mathbf{F}_{\mathrm{coadd}}\delta(x,y).$ (24)

If the coadd process is not a linear combination of the inputs, then
$\mathbf{F}_{\mathrm{coadd}}$ is not a linear operator, and it is no longer
valid to write $I_{\mathrm{coadd}}$ as a convolution of $P_{\mathrm{coadd}}$:

$\displaystyle I_{\mathrm{coadd}}(x,y)$
$\displaystyle=\mathbf{F}_{\mathrm{coadd}}\int
dx^{\prime}dy^{\prime}T(x^{\prime},y^{\prime})\delta(x-x^{\prime},y-y^{\prime})$
$\displaystyle\neq\int
dx^{\prime}dy^{\prime}T(x^{\prime},y^{\prime})\left(\mathbf{F}_{\mathrm{coadd}}\delta(x-x^{\prime},y-y^{\prime})\right)$
$\displaystyle I_{\mathrm{coadd}}(x,y)$ $\displaystyle\neq T(x,y)\otimes
P_{\mathrm{coadd}}(x,y).$ (25)

Figure 4: The effect of using a median coadd on an extended object. The thin
blue, orange, and green curves show a single exponential galaxy profile
convolved by the same three PSF profiles as shown in Figures 1 and 2. The
broad, semi-transparent grey curve is the result of using a median coadd
procedure on these images. The broad, semi-transparent pink curve is the
convolution of the true galaxy profile with the inferred PSF from the median
stellar profile at high signal-to-noise (i.e., the pink profile in Figure 2).
The fact that the grey and pink curves differ shows that even at high signal-
to-noise, nonlinear coaddition algorithms (like median) do not result in a
well-defined PSF.

Thus, when the coadd process is not linear, the coadd image of an extended
source is not the convolution of the true scene by what we would want to call
the PSF. Essentially, the entire concept of what we mean by a “point-spread
function” (i.e., the response of the imaging operator to a delta function, or
“point”) is invalid for a nonlinear imaging operator. Since the coaddition
process is properly included as part of this operator, if that step is
nonlinear, then there is no well-defined PSF. Figure 4 demonstrates this fact
for an exponential galaxy profile using a median coaddition process (using the
same three PSFs as shown in Figures 1 and 2).

##  5 Conclusion

In this paper, we showed a number of choices of coaddition algorithm that
could result in the coadd image not having a well-defined PSF, by which we
mean that there is no PSF function one can define for each location in the
coadd such that the observed image is a convolution of the true scene with
this PSF. In general, for a coadd PSF to be well-defined, the coaddition
procedure must be linear, and the weights assigned to each image in the linear
coadd must be uncorrelated with the true sky scene. For example, both of the
following choices that are sometimes used in real data would result in a coadd
that has a well-defined PSF: (a) a spatially constant weight that is inversely
proportional to the average PSF size in an exposure, and (b) a spatially
varying weight that is inversely proportional to the sky brightness (not
including the astronomical sources).

We note that not having a well-defined PSF is not necessarily a problem. There
are many use cases for coadd images where one does not care whether the PSF is
well-defined – e.g., making color images for showing observations of a
planetary nebula, or possibly detection coadds where one may just want to
maximize total signal-to-noise. But for almost any use case where one will be
making some measurement (e.g., flux or shape) based on the pixel values in the
coadd, it is essential to have a well-defined PSF. In these cases, the
tradeoff between systematic and statistical uncertainties may point towards
choices that are less statistically optimal (e.g., no use of inverse variance
weights that include the object flux) but that reduce systematic errors. Note,
moreover, that this has implications for steps that may often be considered
part of a process that precedes coaddition, such as cosmic ray rejection (see,
for example, the discussion in section 3.3 of Bosch et al. 2018).

We found that the following choices in the coadd algorithm lead to the coadd
failing to have a well-defined PSF:

  * •

All nonlinear coadd algorithms fail to have a well-defined PSF. The coadd
procedure is essentially part of the total imaging process, so if that is not
linear in the inputs, then the final profiles of extended objects are not the
convolution of the true scene with the profile of a delta function ("point"),
thus invalidating the whole idea of a "point-spread function".

  * •

Even for point sources, nonlinear coadd algorithms cause the expectation value
of the observed profile shape to change as a function of signal-to-noise, so
the expected profile is not independent of flux. Thus, they do not have a
well-defined PSF even when only applied to stars.

  * •

Linear coadd algorithms fail to have a well-defined PSF if their weights are a
function of the signal. This means if one is using inverse variance weighting,
one should not include the Poisson noise of the signal as part of that
variance.

  * •

Technically, any weight function that is not spatially constant on each source
image leads to the final image not having a formally well-defined PSF in the
sense of following Eq. (5). In practice this is often not a problem in an
expectation-value sense, so long as the spatial variability of the weights is
uncorrelated with the input signal. In this case, the coadd image follows Eq.
(18). The non-uniform weights can add some additional scatter in the observed
$I_{\mathrm{coadd}}(x,y)$, but not a bias.

Given these constraints, we are left only with linear coadd schemes using
weights that are independent of the signal. If one needs the coadd PSF to be
specifically correct at every location in the image, rather than merely
unbiased, then one should use spatially invariant weights. A strategy that is
usually acceptable in practice is to use weights that vary slowly across the
image, so that the weights are nearly uniform for most objects of interest.
This may be particularly important for science involving two-point functions
of measured quantities, since having a zero expectation value for the bias may
not be sufficient to avoid biasing scientific analyses. Note that we did not
go into detail about matrix-based coadds, but these also qualify as linear, so
they can be perfectly valid, so long as the matrix of weights follows the same
restrictions as outlined above.

## Acknowledgments

We thank Pat Burchat, Johann Cohen-Tanugi, Xiangchong Li, Josh Meyers, Morgan
Schmitz, and Chris Stubbs for useful feedback on early drafts of this paper.
The scripts that produce the plots in this paper are available at
https://github.com/LSSTDESC/DESCNote_CoaddPsf. This paper has undergone
internal review in the LSST Dark Energy Science Collaboration. The internal
reviewers were Gary Bernstein, Johann Cohen-Tanugi, and Ismael Mendoza.

RM and TZ acknowledge the support of the Department of Energy grant DE-
SC0010118. MJ is partially supported by Department of Energy grant DE-
SC0007901.

Author contributions to this work are as follows: RM synthesized telecon
discussions and material from the literature to produce initial draft Note,
discussed with co-authors, edited. MJ srote significant portions of the text,
including 2.1, part of 3.2, 4.3, figure 4, plus edits and word smithing
throughout. RHL contributed to the study which led to this paper, and provided
comments/suggestions during the writing of the paper. JB did some of the early
work (on the HSC pipelines) that first identified the problems more fully
described in this paper, and participated in discussions and editing of early
drafts. AK contributed to some of clarification and general feedback to the
publication, and created a notebook to numerically verify some of the key
equations (Eqs. 10-11). MDM created the median coadd bias quantification in
Figure 3, as well as supporting text in Section 4.1. TZ double-checked the
nonlinear weighting coadd with more specifications, and provided feedback on
the impact of median coadd results.

The DESC acknowledges ongoing support from the Institut National de Physique
Nucléaire et de Physique des Particules in France; the Science & Technology
Facilities Council in the United Kingdom; and the Department of Energy, the
National Science Foundation, and the LSST Corporation in the United States.
DESC uses resources of the IN2P3 Computing Center (CC-IN2P3–Lyon/Villeurbanne
- France) funded by the Centre National de la Recherche Scientifique; the
National Energy Research Scientific Computing Center, a DOE Office of Science
User Facility supported by the Office of Science of the U.S. Department of
Energy under Contract No. DE-AC02-05CH11231; STFC DiRAC HPC Facilities, funded
by UK BIS National E-infrastructure capital grants; and the UK particle
physics grid, supported by the GridPP Collaboration. This work was performed
in part under DOE Contract DE-AC02-76SF00515.

## References

  * Annis et al. (2014) Annis, J., Soares-Santos, M., Strauss, M. A., et al. 2014, ApJ, 794, 120
  * Bosch et al. (2018) Bosch, J., Armstrong, R., Bickerton, S., et al. 2018, PASJ, 70, S5
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111
  * Lauer (1999) Lauer, T. R. 1999, Publications of the Astronomical Society of the Pacific, 111, 227 – 237, cited by: 82
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration, Abell, P. A., Allison, J., et al. 2009, arXiv e-prints, arXiv:0912.0201
  * Mandelbaum (2018) Mandelbaum, R. 2018, ARA&A, 56, 393
  * Rowe et al. (2011) Rowe, B., Hirata, C., & Rhodes, J. 2011, ApJ, 741, 46
  * Rowe et al. (2015) Rowe, B. T. P., Jarvis, M., Mandelbaum, R., et al. 2015, Astronomy and Computing, 10, 121
  * Zackay & Ofek (2017) Zackay, B., & Ofek, E. O. 2017, ApJ, 836, 187