# Cycle Text-to-Image GAN with BERT

Trevor Tsue
ttsue
Computer Science Dept. &Jason Li
jasonkli
Computer Science Dept. &Samir Sen
samirsen
Computer Science Dept

(9 June 2019)

###### Abstract

We explore novel approaches to the task of image generation from their
respective captions, building on state-of-the-art GAN architectures.
Particularly, we baseline our models with the Attention-based GANs that learn
attention mappings from words to image features. To better capture the
features of the descriptions, we then built a novel cyclic design that learns
an inverse function to maps the image back to original caption. Additionally,
we incorporated recently developed BERT pretrained word embeddings as our
initial text featurizer and observe a noticeable improvement in qualitative
and quantitative performance compared to the Attention GAN baseline. 111Cycle
Image GAN Github

##  1 Introduction

(a) Text to Image

(b) Text to Image to Text

Figure 1: Text to Image Generation Tasks

The goal of the text-to-image task is to generate realistic images given a
text description. This problem has many possible applications ranging from
computer-aided design to art generation Xu et al. (2017). Moreover, this
multimodal problem is an interesting and important task in natural language
understanding because it connects language to an understanding of the visual
world.

The problem can be naturally decomposed into two parts: embedding the text
into a feature representation that captures relevant visual information and
using that representation to generate an realistic image that corresponds to
the text. This problem is particularly challenging for several reasons. For
one, there is the issue of a domain gap between the text and image feature
representations. Furthermore, there are a myriad of legitimate images that
correspond to a single text description and part of the goal is to be able to
capture the diversity in plausible images Reed et al. (2016).

Recent advances in the field of deep learning have made significant strides in
this challenging task. In particular, recurrent architectures, such as LSTMs,
can be used to learn feature representations from text and generative
adversarial networks (GANs) can be used to create images conditioned on
information. Nonetheless, the aforementioned challenges leave the text-to-
image task an open problem. Adding to those problems, GANs, despite having
widespread success in generative learning, often produce lower-resolution
images, Zhang et al. (2016), lack diversity in image generation, and fail to
capture intricate details.

The goal of our work is to explore and compare state-of-the-art methods for
addressing some of these problems. These ideas include: stacking multiple GANs
to sequentially learn higher resolution images; applying attention mechanisms
to “focus” the generator on important parts of the text and to help bridge the
domain gap; and unifying text-to-image and image-to-text in a single model. In
addition to evaluating these ideas for diversity and quality of generated
images, we also investigate the effect of using a pre-trained language model
for contextualized word embeddings. Pretrained language models, such as BERT
Devlin et al. (2018) and ELMO Peters et al. (2018), have revolutionized NLP as
an effective means of transfer learning, similar to the impact of ImageNet on
the field of computer vision. As such, we sought to explore the potential
benefit of these pretrained embeddings since most current approaches learn the
word embeddings from scratch.

##  2 Related Work

Generative adversarial networks (GANs) are the most widely used model in
generative learning. Originally proposed to generate realistic images, GANs
consist of two neural networks, a generator and a discriminator, that
effectively compete with one another in a zero-sum game. The discriminator
attempts to distinguish real and fake images, while the generator tries to
create images that fool the discriminator into classifying them as real
Goodfellow et al. (2014).

Numerous works have built off of the basic idea. Some of these idea include
the Conditional GANs, which pass a class label to both the generator and
discriminator, unlike the original GAN where the generator creates an image
solely from noise Mirza and Osindero (2014). GANs have also been used for
style transfer between image domains. In this formulation, the generator is
passed an image from a source domain and tries to fool the discriminator into
thinking it is from the target domain. An extension to this idea is the
CycleGAN, which learns to transfer from source to target and target back to
source to ensure consistency and to stabilize training. This setup also
ensures that important latent features are captured so that the source image
can be reconstructed from the generated one Zhu et al. (2017). As we will see
next, these ideas have a natural extension to the text-to-image problem.

Reed et al. describe the first fully differentiable, end-to-end model that
learns to construct images from text, building on conditional GANs Reed et al.
(2016). They use a character-level convolutional-recurrent network to encode
the input text. A fully-connected (FC) layer with a Leaky ReLU activation
embeds the encoding into a lower dimension before concatenation with input
noise drawn from a standard Gaussian, which is then fed into the generator.
Instead of just training on pairs of either real images/matching text or fake
images/matching text, they also train using pairs of real images with
mismatched text. This is to encourage the discriminator to not only generate
realistic images regardless of the text, but also to create realistic images
that match the text.

Zhang et al. build upon this work in “StackGAN: Text to Photo-realistic Image
Synthesis with Stacked Generative Adversarial Networks”. Inspired by other
works that use multiple GANs for tasks such as scene generation, the authors
used two stacked GANs for the text-to-image task Zhang et al. (2016). The
motivating intuition is that the Stage-I GAN produces a low-resolution outline
of the desired image and the Stage-II GAN fills in the details of the sketch.
In addition to the stacked architecture, they also propose a novel
augmentation technique to address the aforementioned interpolation issues.
Rather than using the text embedding directly (following a single FC layer) as
done in Reed et al., they instead use an FC layer to produce a mean and
variance before sampling from the normal distribution. This sampling augments
the data and increases the robustness.

Another architecture developed by Xu et al. drew on the widespread success of
attention-based models, particularly in NLP for tasks such as machine
translation and image captioning. Their model, termed AttnGAN, introduces
several novel ideas to the text-to-image task related to attention Xu et al.
(2017). Unlike previous approaches that focus on sentence-level encodings,
AttnGAN extracts both sentence-level and word-level features from a
bidirectional LSTM. In the stacked generator stages, multiplicative attention
is performed over the encoded word vectors so the model can learn which words
to attend to at each step. Finally, they add a Deep Attentional Multimodal
Similarity Model, which is constructed to learn an attention-based matching
score between the image-sentence pairs.

Finally, Gorti et al. incorporate the ideas of stacking, attention, and cycle
consistency in their state-of-the-art model, MirrorGAN Qiao et al. (2019).
Influenced by the CycleGAN architecture, the model adds an image-to-text
component which acts as a sanity check that the image generated is indeed
semantically consistent with the input caption text. The results demonstrate
MirrorGAN’s ability to train networks that can generate both higher quality
images as well as image details which are semantically consistent with a
provided caption (and in comparison with the true image for a given caption).

MirrorGAN is the culmination of the work on the text-to-image problem.
Nonetheless, there is room for improvement. The works up until this point use
word embeddings trained from scratch. With the advent of pretrained language
models such as ELMO or BERT, a possible extension is to initialize the
embeddings with deep, contextualized word vectors derived from BERT or ELMO.
In this work, we explore the effect of using BERT-derived word vectors.

##  3 Data

We used the 2011 Caltech-UCSD Birds 200 dataset (CUB-200), which contains
11,788 images of 200 different types of birds and is a widely used benchmark
for text-to-image generation Wah et al. (2011). These images provide a
boundary box and vary in size. Additionally, we have 10 text descriptions of
the dataset downloaded from a github repository that serve as the text
descriptions of the generated images 222taoxugit AttnGAN .

##  4 Methods

###  4.1 Data Preprocessing

We preprocess this data according the precedence set by StackGAN++ Zhang et
al. (2017). This includes cropping all images to ensure all bounding boxes
have at least a 0.75 object-image size ratio and then downsampled to 64x64,
128x128, 256x256. Then, the data is split into class disjoint train and test
sets.

###  4.2 Models

####  4.2.1 AttnGAN

Figure 2: Attention GAN Architecture

This first model combines both elements of the stack GAN Zhang et al. (2016)
and attention Xu et al. (2017) . This attention GAN first embeds the caption
and runs them through a LSTM, generating both word and sentence vectors. Using
the Conditioning Augmentation first proposed in the StackGANSalimans et al.
(2016), we create a mean and variance from the sentence embedding via a fully-
connected layer. We use this mean and variance to parameterize a normal
distribution from which a sentence embedding sample is generated to pass into
the GAN. This is used for regularization and to promote manifold smoothness.
Additionally, we concatenate Gaussian noise to this new sentence embedding
sample and pass into the generator.

With the StackGAN architecture, we stack three generators together, generating
64x64, 128x128, and 256x256, respectively. Additionally, for the second and
third generator, we pass the image and the word embeddings through an
attention module to pass into the next generator Xu et al. (2017). Each of
these generators has a corresponding discriminator that take in both the
original sentence embedding and the image. Finally, the 256x256 image is
passed through an image encoder to generate local image features (a 17x17
feature map). These image features from the image encoder and word features
from the text encoder combine to form the Deep Attentional Multimodal
Similarity Model (DAMSM) and trained with an attention loss Xu et al. (2017).
For stability, we pretrained this DAMSM model.

####  4.2.2 CycleGAN

Figure 3: Cycle GAN Architecture

Our CycleGAN combines the attention GAN and the original cycle GAN approach
Zhu et al. (2017). By adding an RNN conditioned on the image features and the
embedded captions, we attempt to return to text with the Semantic Text
REgeneration and Alignment Module (STREAM) Qiao et al. (2019). By learning
this transition to the original text domain, we allow our images to better
represent our captions, as they must hold the latent information to recreate
the original caption Qiao et al. (2019) 333komiya-m MirrorGAN. Additionally,
we added the pretrained BERT encoding transformers, which we use instead of
the standard word embeddings.

####  4.2.3 BERT

Pretrained word vectors are a common component of many of NLP models. However,
until recently, one primary limitation of these word vectors was that they
only allowed for one context-independent embedding. One of the biggest game-
changers in recent NLP research is the advent of deep, contextualized word
vectors. These vectors are derived from the internal states of deep,
pretrained language models trained on massive corpuses of text and use the
entire sequence to embed each word, not just the word itself. A key premise to
this idea was previous research that showed different layers of an LSTM
language model captured different information, such as part of speech at the
lower levels and context at the higher layers Peters et al. (2018)

Peters et al. were the first to introduce this idea with their language model,
ELMO.. Their model was a deep, bidirectional LSTM with character-level
convolutions. This pretrained model could then be used for more specific
tasks, where each word vector is computed as the (learned) weighted sum of the
hidden states of the LSTM, using the entire input sequence as input. Then,
tasks like sentiment classification could be done by simply adding a fully
connected layer on top of ELMO Peters et al. (2018). Devlin et al. expanded on
this work with the BERT model that replaces the bidirectional LSTM with a
bidirectional Transformer Devlin et al. (2018). The Transformer is another
recent innovation in NLP that replaces the recurrent nature of LSTMs with
positional encoding and blocks of self-attention, layer normalization, and
fully connected layers Vaswani et al. (2017). This architecture has become the
de facto model for NLP, replacing LSTMs and standard RNNs in many cases. The
BERT model has been widely used to achieve state-of-the-art results in
challenging tasks such as Question-Answer (QA). In this work, we used a
pretrained BERT model to obtain our embeddings and pass it through a fully
connected layer before continuing in the CycleGAN architecture.

###  4.3 Loss

For a generator $G_{i}$ and the corresponding discriminator $D_{i}$, we have
the following loss function which combines both a conditional and
unconditional loss (conditioned on the sentence embedding):

$\displaystyle\mathcal{L}_{G_{i}}=-\frac{1}{2}\mathbb{E}_{\hat{x}_{i}\sim
p_{G_{i}}}[\log(D_{i}(\hat{x}_{i}))]$
$\displaystyle-\frac{1}{2}\mathbb{E}_{\hat{x}_{i}\sim
p_{G_{i}}}[\log(D_{i}(\hat{x}_{i},\bar{e}))]$

where $\hat{x}_{i}$ is the generated image and $\bar{e}$ is the sentence
embedding.

Then, we have the following discriminator loss:

$\displaystyle\mathcal{L}_{D_{i}}=-\frac{1}{2}\mathbb{E}_{x_{i}\sim
p_{data_{i}}}[\log(D_{i}(x_{i}))]$
$\displaystyle-\frac{1}{2}\mathbb{E}_{\hat{x}_{i}\sim
p_{G_{i}}}[\log(1-D_{i}(\hat{x}_{i}))]$
$\displaystyle-\frac{1}{2}\mathbb{E}_{x_{i}\sim
p_{data_{i}}}[\log(D_{i}(x_{i}))]$
$\displaystyle-\frac{1}{2}\mathbb{E}_{\hat{x}_{i}\sim
p_{G_{i},\bar{e}}}[\log(1-D_{i}(\hat{x}_{i},\bar{e}))]$

####  4.3.1 AttnGAN Loss

With the word embeddings matrix $e$ and the image embeddings $v$, we calculate
a similarity score between the sentence:

$\displaystyle s=e^{\intercal}v$ (1)

We create a $s\in\mathbb{R}^{T\times 289}$ with $T$ as the number of words in
the sentence and 289 referring to a flattened version of the 17x17 image
feature map. We then normalize the similarity matrix

$\displaystyle\bar{s}_{ij}=\frac{\exp(s_{ij})}{\sum_{k=1}^{T}\exp(s_{kj})}$
(2)

We build a context vector $c_{i}$, where that represent the image regions that
relate to the $i$th word in the sentence.

$\displaystyle c_{i}=\sum_{j=1}^{289}\alpha_{j}v_{j}\text{, where
}\alpha_{j}=\frac{\exp(\gamma\bar{s}_{ij})}{\sum_{k=1}^{289}\exp(\bar{s}_{ik})}$
(3)

where $\gamma$ is a hyperparameter to pay attention to certain features in the
regions. We then have an attention-driven image-text matching score matching
the entire image $Q$ to the whole text description $D$ that utilizes the
cosine similarity
cosine$(c_{i},e_{i})=\frac{c_{i}^{\intercal}e_{i}}{\|c_{i}\|\|e_{i}\|}$

$\displaystyle
R(Q,D)=\log\Big{(}\sum_{i=1}^{T}\exp(\gamma\text{cosine}(c_{i},e_{i}))\Big{)}$
(4)

We have the DAMSM probability between the different image-sentence pairs in
the batch

$\displaystyle P(D_{i}|Q_{i})=\frac{\exp(\gamma
R(Q_{i},D_{i}))}{\sum_{j=1}^{M}\exp(\gamma R(Q_{i}D_{j}))}$ (5) $\displaystyle
P(D_{i}|Q_{i})=\frac{\exp(\gamma R(Q_{i},D_{i}))}{\sum_{j=1}^{M}\exp(\gamma
R(Q_{j}D_{i}))}$ (6)

Therefore, our DAMSM combines the following

$\displaystyle\mathcal{L}_{1}^{w}=-\sum_{i=1}^{M}\log P(D_{i}|Q_{i})$ (7)
$\displaystyle\mathcal{L}_{2}^{w}=-\sum_{i=1}^{M}\log P(Q_{i}|D_{i})$ (8)

We also define $\mathcal{L}_{1}^{s}$ and $\mathcal{L}_{2}^{s}$ the same as
above but instead substituting $\bar{e}$ for $e$. Combining everything, we
have the final loss of the attention generator

$\displaystyle\mathcal{L}_{DAMSM}=\mathcal{L}_{1}^{w}+\mathcal{L}_{2}^{w}+\mathcal{L}_{1}^{s}+\mathcal{L}_{2}^{s}$
(9)
$\displaystyle\mathcal{L}=\mathcal{L}_{G}+\lambda\mathcal{L}_{DAMSM}\text{,
where }\mathcal{L}_{G}=\sum_{i=1}^{3}\mathcal{L}_{G_{i}}$ (10)

####  4.3.2 CycleGAN Loss

In addition to the loss of the AttnGAN, we add an additional cross entropy
loss to correctly predict the output word in the caption recreation

$\displaystyle\mathcal{L}_{CE}=-\frac{1}{M}\sum_{i=1}^{M}\sum_{c=1}^{|V|}y_{c}^{(i)}\log(\hat{y}_{c}^{(i)})$
$\displaystyle\mathcal{L}=\mathcal{L}_{G}+\lambda\mathcal{L}_{DAMSM}+\lambda\mathcal{L}_{CE}$

where $M$ represents the batch size, $|V|$ is the size of the vocab,
$y_{c}^{(i)}$ is the binary label of the $c$-th class of the $i$-th example,
and $\hat{y}_{c}^{(i)}$ is the model probability output of the $c$-th class of
the $i$-th example.

###  4.4 Evaluation Metrics

####  4.4.1 Inception Score

To assess our models, we use the Inception score, which is a widely used ad
hoc metric for generative models Salimans et al. (2016). The Inception score
uses a pretrained Inception model that is fine-tuned to the specific dataset
being used. The Inception score is computed by exponentiating the KL-
divergence between the conditional distribution p(y $|$ x) and marginal
distribution p(y), where y is the class label predicted by the Inception model
and x is a generated sample. The intuition is that a good generative model
should produce images with a conditional label distribution that has low
entropy relative to the marginal distribution. In other words, we want images
that can be easily classified into a category by the model but also create
images that belong to many different classes.

$\displaystyle
D_{KL}(P\|Q)=-\sum_{x\in\mathcal{X}}P(x)\log\Big{(}\frac{Q(x)}{P(x)}\Big{)}$
$\displaystyle IS(G)=\exp\Big{(}\mathbb{E}_{\mathbb{x}\sim
p_{G}}D_{KL}(p(y|\mathbb{x})\|p(y))\Big{)}$

The score rewards images that have greater variety and has been shown to be
well-correlated with human evaluations of realistic quality. We randomly
select 20 captions for each class and use our trained model to generate
images, which is then fed into the Inception model to generate the
distributions and to compute the score.

####  4.4.2 Mean Opinion Score (MOS)

Nonetheless, the Inception score cannot capture how well the generated images
reflect accurate conditioning on the input text. Thus, we have humans examine
the perceptual quality of images as well as their correspondence to the input
task with the Mean Opinion Score Ledig et al. (2016). Specifically, we asked
$n=10$ subjects to rate the quality of images on a scale from 1 (poor quality)
to 5 (high quality). We showed them 20 images from the ground truth, 20 from
the AttnGAN, and 20 from the CycleGAN, along with corresponding captions, in
random order, and averaged them to report the MOS.

##  5 Results

Figure 4: Generated Images from the models

We trained the AttnGAN model over 100 epochs using Adam optimization to train
the generator and all the discriminators. CycleGAN was trained over 100 epochs
using the same generator and discriminator optimizers with betas of 0.5 and
0.999. For the AttnGAN, we pretrained the DAMSM architecture for 200 epochs.
For the CycleGAN, we pretrained the STREAM architecture for 100 epochs. We
used pretrained BERT embeddings only for the CycleGAN implementation, while
initialized randomly initialized embeddings which were trained in the AttnGAN.

We report both the inception v3 scores, computed as a average measure of
divergence with the true distribution of bird images of the generated test
outputs for each of the models as well as the qualitative MOS scores from peer
judges, reported below. We see that the CycleGAN trained with BERT embeddings
had the strongest performance overall across the proposed metrics, and display
generated samples from our model along with their representative ground truth
image labels.

Inception Score Model Epoch 0 Epoch 50 Epoch 100 Ground Truth 11.63 - -
AttnGAN 0.94 2.78 3.92 CycleGAN w/ BERT 1.05 5.48 5.92

Figure 5: Inception Scores of Models

Model MOS (n=10) Ground Truth 4.7 AttnGAN 3.6 CycleGAN w/ BERT 3.9

Figure 6: Mean Opinion Score (MOS) of Models with $n=10$ subjects

We save model weights for the CycleGAN model with BERT text features every 25
epochs and compute Inception scores on the validation set during training. In
Figure 7, we observe the CycleGAN inception score leveling, but still
increasing as we approach 100 epochs of training.

Figure 7: CycleGAN Inception Scores

##  6 Discussion

Examining several images output from the AttnGAN and CycleGAN (with BERT) in
Figure 4, we can see some clear improvements from AttnGAN to CycleGAN. For
one, the CycleGAN model generally produces clearer, more realistic looking
images relative to AttnGAN. Furthermore, the CycleGAN model appears to be more
precise with respect to details. We can see that for the AttnGAN model, the
colors are occasionally incorrect (presence of red in the top image and brown
instead of grey in the bottom). Additionally, AttnGAN images lack the level of
detail in the beak present in the CycleGAN model.

We found the CycleGAN with pretrained BERT embeddings was able to outperform
AttnGAN on the test set in both the Inception score (Figure 5) and in the Mean
Opinion Score (Figure 6). In particular, with respect to Inception scores, we
found that CycleGAN with BERT was able to reach a higher score and train
significantly faster, as indicated by the scores at 50 epochs (5.48 vs 2.78)
and at 100 epochs (5.92 vs 3.92). In general, a higher Inception score
reflects greater variety as well as distinctly capturing unique features, but
we note that an ideal quantitative metric remains elusive for this task,
particularly in capturing the correspondence between the image and caption.
The higher performance of the CycleGAN with BERT on qualitative, human
evaluation indicates some level of improvement in the image-text
correspondence. One limitation of our work is that we were not able to train
until convergence of the Inception score for comparison in the limit, which we
note as a possible avenue for future work.

##  7 Conclusion

In this paper, we investigate the text-to-image generation task by
experimenting with state-of-the-art architectures and incorporating the latest
innovations in NLP, namely, the use of deep contextualized word vectors from
pretrained languagew models, such as BERT. Our baseline model is the AttnGAN,
which utilizes several key features, including stacking of GANs to
progressively learn more detail at higher resolution and attention over word
features, a technique that has found widespread success in a variety of NLP
tasks. Our main model adds two additional features: a cyclic architectures
that adds the image-to-text task in addition to the test-to-image task, and
the use of contexualized word embeddings from a pretrained BERT. Through both
qualitative and quantitative metrics, we found that the addition of these
features showed improved generation of images conditioned on the text and had
faster learning.

For future work, it would be useful to train the models longer until
convergence in some metric (such as Inception score) is reached for complete
analysis. In addition, we did not do any hyperparameter tuning due to time
constraints and thus further improvement may be found through a hyperparameter
search. Further, with more time, it would also be interesting to perform
ablation studies on our full model to show the additional gain, if any,
achieved from adding only BERT or only the cyclic architecture to AttnGAN.

##  8 Acknowledgements

We would like to thank the instructors and TAs for designing and running the
course. In particular we would like to thank Ignacio Cases for guiding us on
this project.

##  9 Authorship Statement

Trevor Tsue: Coded the AttnGAN and CycleGAN, trained models, made architecture
diagrams, wrote architecture and loss.

Jason Li: Performed literature review, proposed and coded/integrated BERT
extension with CycleGAN; wrote intro, related work, discussion, conclusion

Samir Sen: AttnGAN implementation and training. Worked on incorporating BERT
encodings within the attention GAN text featurization and built inception
network for capturing key metrics across models. Data cleaning and
visualization. Results, discussion, abstract.

## References

  * Devlin et al. (2018) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2018.  BERT: pre-training of deep bidirectional transformers for language understanding.  _CoRR_ , abs/1810.04805.
  * Goodfellow et al. (2014) Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio. 2014.  Generative adversarial nets.  In Z. Ghahramani, M. Welling, C. Cortes, N. D. Lawrence, and K. Q. Weinberger, editors, _Advances in Neural Information Processing Systems 27_ , pages 2672–2680. Curran Associates, Inc.
  * Ledig et al. (2016) Christian Ledig, Lucas Theis, Ferenc Huszar, Jose Caballero, Andrew P. Aitken, Alykhan Tejani, Johannes Totz, Zehan Wang, and Wenzhe Shi. 2016.  Photo-realistic single image super-resolution using a generative adversarial network.  _CoRR_ , abs/1609.04802.
  * Mirza and Osindero (2014) Mehdi Mirza and Simon Osindero. 2014.  Conditional generative adversarial nets.  _CoRR_ , abs/1411.1784.
  * Peters et al. (2018) Matthew E. Peters, Mark Neumann, Mohit Iyyer, Matt Gardner, Christopher Clark, Kenton Lee, and Luke Zettlemoyer. 2018.  Deep contextualized word representations.  _CoRR_ , abs/1802.05365.
  * Qiao et al. (2019) Tingting Qiao, Jing Zhang, Duanqing Xu, and Dacheng Tao. 2019.  Mirrorgan: Learning text-to-image generation by redescription.  _CoRR_ , abs/1903.05854.
  * Reed et al. (2016) Scott E. Reed, Zeynep Akata, Xinchen Yan, Lajanugen Logeswaran, Bernt Schiele, and Honglak Lee. 2016.  Generative adversarial text to image synthesis.  _CoRR_ , abs/1605.05396.
  * Salimans et al. (2016) Tim Salimans, Ian J. Goodfellow, Wojciech Zaremba, Vicki Cheung, Alec Radford, and Xi Chen. 2016.  Improved techniques for training gans.  _CoRR_ , abs/1606.03498.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin. 2017.  Attention is all you need.  _CoRR_ , abs/1706.03762.
  * Wah et al. (2011) C. Wah, S. Branson, P. Welinder, P. Perona, and S. Belongie. 2011.  The Caltech-UCSD Birds-200-2011 Dataset.  Technical Report CNS-TR-2011-001, California Institute of Technology.
  * Xu et al. (2017) Tao Xu, Pengchuan Zhang, Qiuyuan Huang, Han Zhang, Zhe Gan, Xiaolei Huang, and Xiaodong He. 2017.  Attngan: Fine-grained text to image generation with attentional generative adversarial networks.  _CoRR_ , abs/1711.10485.
  * Zhang et al. (2016) Han Zhang, Tao Xu, Hongsheng Li, Shaoting Zhang, Xiaolei Huang, Xiaogang Wang, and Dimitris N. Metaxas. 2016.  Stackgan: Text to photo-realistic image synthesis with stacked generative adversarial networks.  _CoRR_ , abs/1612.03242.
  * Zhang et al. (2017) Han Zhang, Tao Xu, Hongsheng Li, Shaoting Zhang, Xiaogang Wang, Xiaolei Huang, and Dimitris N. Metaxas. 2017.  Stackgan++: Realistic image synthesis with stacked generative adversarial networks.  _CoRR_ , abs/1710.10916.
  * Zhu et al. (2017) Jun-Yan Zhu, Taesung Park, Phillip Isola, and Alexei A. Efros. 2017.  Unpaired image-to-image translation using cycle-consistent adversarial networks.  _CoRR_ , abs/1703.10593.


# A general method to reconstruct strong gravitational lenses based on the
singular perturbative approach.

Alard, C.,
IAP, 98bis Boulevard Arago, Paris

###### Abstract

The number of gravitational arcs systems detected is increasing quickly and
should even increase at a faster rate in the near future. This wealth of new
gravitational arcs requires the development of a purely automated method to
reconstruct the lens and source. A general reconstruction method based on the
singular perturbative approach is proposed in this paper. This method
generates a lens and source reconstruction directly from the gravitational arc
image. The method is fully automated and works in two steps. The first step is
to generate a guess solution based on the circular solution in the singular
perturbative approach. The second step is to break the sign degeneracy and to
refine the solution by using a general source model. The refinement of the
solution is conducted step by step to avoid the source-lens degeneracy issue.
One important asset of this automated method is that the lens solution is
written in universal terms which allows the computation of statistics.
Considering the large number of lenses which should be available in the near
future this ability to compute un-biased statistics is an important asset.

###### keywords:

gravitational lensing: strong - (cosmology:) dark matter

††pubyear: 2023††pagerange: A general method to reconstruct strong
gravitational lenses based on the singular perturbative
approach.–LABEL:lastpage

##  1 Introduction

There are basically two main issues with the automated reconstruction of
gravitational arcs: (1) what model to choose for a given gravitational arc ?
(2) Given the large parameter space for the models, how to identify some first
guess for the solution ? We will show in this paper that the singular
perturbative approach offers an efficient solution to these two issues. The
paper starts with an introduction to the perturbative singular approach, and
continues with a discussion on the proper referential for the reconstruction
of the lens. In this referential a first guess for the solution is constructed
based on the circular source model in the perturbative approach (see Alard
(2007) Eq. (12)). Once this first guess is reconstructed the next step will be
to use a general chi-square minimization method to refine the former guess and
reach the optimal $\chi_{2}$ minimum. The degeneracy of the first guess is
naturally broken by constructing a refined non circularly symmetric source
model. The reconstruction is completed by increasing the potential Fourier
order to the required minimal order. The first part of the paper will recall
the basics of the singular perturbative approach. This will be followed by a
description of the methods and associated equations. The presentation will be
illustrated with a number of practical applications by reconstructing
simulated gravitational arcs images.

##  2 Basics of the singular perturbative approach in gravitational lensing.

We define the local coordinates ${\bf r}$ and source coordinates ${\bf
r_{S}}$, in a lens potential $\phi$ the lens equation reads,

${\bf r_{S}}={\bf r}-\nabla\phi$ (1)

The basic idea of the singular perturbative approach is to consider that the
un-perturbed situation corresponds to the image of a point at the center of a
circularly symmetric potential. Due to its symmetry the image of a point in a
circularly symmetric potential is a circle corresponding to the Einstein
circle ($R_{E}$). For convenience we re-normalize the coordinate system to
obtain $R_{E}=1$. As consequence the unperturbed situation is an infinity of
points located on the Einstein circle. The choice of this un-perturbed
situation allow us to deal with the high non-linearity of strong gravitational
lenses in the angular direction. In effect at any angular position $\theta$ an
un-perturbed point is present. On the other hand, the perturbation in the
radial direction is generally small for gravitational arcs, which allow us to
work in the vicinity of the Einstein circle. It is also assumed that strong
gravitational arcs are due to small perturbation of a circular
potential.Consequently in this framework we have,

$\begin{cases}\phi(r,\theta)=\phi_{0}(r)+\epsilon\psi(r,\theta)\\\
r=1+\epsilon dr\\\ {\bf r_{S}}=\epsilon{\bf r_{s}}\end{cases}$ (2)

In Eq. (2) $\epsilon$ is a small number. By introducing Eq. (2) in the lens
equation Eq. (1) the perturbative lens equation is obtained (see Alard (2007)
Eq. (8) and (7)),

$\begin{cases}\mathbb{r}_{s}=(\kappa_{2}dr-
f_{1})\mathbb{u}_{r}-\frac{df_{0}}{d\theta}\mathbb{u}_{\theta}\\\
f_{n}(\theta)=\frac{1}{n!}\left[\frac{\partial^{n}\psi}{\partial
r^{n}}\right]_{r=1}\\\
\kappa_{2}=1-\left[\frac{d^{2}\phi_{0}}{dr^{2}}\right]_{r=1}\end{cases}$ (3)

It is useful to include the source impact parameter ${\bf r_{0}}$ in the
definition of the perturbative fields $f_{1}$ and $\frac{df_{0}}{d\theta}$ by
defining the new fields $\tilde{f}_{1}$ and $\frac{d\tilde{f}_{0}}{d\theta}$,

$\tilde{f}_{i}=f_{i}+{\bf r_{0}u_{r}}~{}~{}~{}~{}{\rm where}~{}i=0,1$ (4)

In the continuation of this paper we will assume the impact parameters are
included in the field and will not use the ${\tilde{f}_{i}}$ symbol to lighten
the notation and will use directly the $f_{i}$ notation.

###  2.1 The solution for a circular source in the perturbative approach.

An interesting application is to consider the image of a round source contour
with radius $R_{0}$ in the singular perturbative theory. This model is simple
with a direct analytical solution for the images. In effect in this case
solving Eq. (3) is straightforward and reduces to a second order equation (see
Alard (2007), Eq. 12).

$dr=\frac{1}{\kappa_{2}}\left[f_{1}\pm\sqrt{R_{0}^{2}-\left(\frac{df_{0}}{d\theta}\right)^{2}}\right]$
(5)

It is interesting to note that in this circular source model the two
perturbative fields have simple physical meanings. The field $f_{1}$
represents the mean position of the image for each angular position, while the
field $\frac{df_{0}}{d\theta}$ is directly related to the image width.

###  2.2 Relation to the multipole expansion

A final and interesting point about the singular perturbative expansion is
that the perturbative fields are related to the multipole expansion on the
Einstein circle $R_{E}$ (here re-scale at $R_{E}=1$). The multipole expansion
of the potential at $R_{E}=1$ reads:

$\psi=-\sum_{n}{\frac{a_{n}(r)}{r^{n}}\cos n\theta+\frac{b_{n}(r)}{r^{n}}\sin
n\theta+c_{n}(r)r^{n}\cos n\theta+d_{n}(r)r^{n}\sin n\theta}$ (6)

The coefficients, $a_{n}(r),b_{n}(r),c_{n}(r),d_{n}(r)$ of the potential
multipole expansion (Eq. 6) are related to the surface density $\Sigma$ by the
following equations,

$\begin{cases}a_{n}=\frac{1}{2\pi
n}\int_{0}^{\pi}\int_{0}^{r=1}\Sigma(u,v)\cos nv\ u^{n+1}dudv\\\
b_{n}=\frac{1}{2\pi n}\int_{0}^{\pi}\int_{0}^{r=1}\Sigma(u,v)\sin nv\
u^{n+1}dudv\\\ c_{n}=\frac{1}{2\pi
n}\int_{0}^{\pi}\int_{r=1}^{\infty}\Sigma(u,v)\cos nv\ u^{1-n}dudv\\\
d_{n}=\frac{1}{2\pi n}\int_{0}^{\pi}\int_{r=1}^{\infty}\Sigma(u,v)\sin nv\
u^{1-n}dudv\\\ \end{cases}$ (7)

This multipole expansion of Eq’s (6) and (LABEL:Eq_pot_multi2) is related to
the singular perturbative expansion by Eq. (8) (see Alard (2009) for the
details of the calculations.)

$\begin{cases}f_{1}=\sum_{n}{n(a_{n}-c_{n})\cos n\theta+n(b_{n}-d_{n})\sin
n\theta}\\\ \frac{df_{0}}{d\theta}=\sum_{n}{-n(b_{n}+d_{n})\cos
n\theta+n(a_{n}+c_{n})\sin n\theta}\end{cases}$ (8)

It is interesting to note that Eq. (8) implies a direct relation between the
Fourier expansion of the perturbative fields and the multipole expansion of
the potential on the Einstein circle. This also implies that in this theory
the inner and outer contribution to the lens potential can be separated. An
example of a reconstruction with a demonstrative separation of the inner and
outer potential contribution can be found in Alard (2010)).

##  3 The method

###  3.1 Estimating the center of the perturbative reconstruction.

A basic problem in the perturbative reconstruction of the lens is that the
actual center of the unperturbed circular distribution is unknown. In most
case the lens itself can be observed and one could use the center of the light
distribution in the lens as the center for the reconstruction. However even if
this assumption seems reasonable it might be wrong in some cases. It is more
interesting to develop a general method without making assumptions and find
the center as a parameter of the reconstruction. We will start by considering
that we are working in a referential which is slightly shifted form the
optimal (true) referential. The effect of the small shift between the two
referential on the perturbative fields will be estimated using a perturbative
expansion.

####  3.1.1 Effect of a small shift on the perturbative fields.

Let’s define the offset from the optimally centered referential,
${\bf\delta_{0}}$. We will assume that, $|\bf\delta_{0}|\ll 1$ or
equivalently, that $|\bf\delta_{0}|$ is of order $\epsilon$ with $\epsilon\ll
1$,

$\bf\delta_{0}=\epsilon\bf\delta$ (9)

The coordinates in the optimal referential and offset referential are
respectively $\bf r$ and $\bf u$. The two coordinates are related by,

${\bf r=u+\delta_{0}}$ (10)

The lens equation (Eq. (1)) is transformed using Eq. (10) leading to,

${\bf r_{S}=u+\delta_{0}-\nabla\phi\left(u+\delta_{0}\right)}$ (11)

By expanding Eq. (11) to first order in $\epsilon$ with the help of Eq. (23)
we obtain,

${\bf r_{S}=u+\delta_{0}-\nabla\left(\phi(u)+\delta_{0}.\nabla\phi(u)\right)}$
(12)

Considering that $\phi=\phi_{0}+\epsilon\psi$ (see Eq. 2) and again expanding
to first order in $\epsilon$ we obtain,

$\begin{cases}{\bf r_{S}=u-\nabla\left(\tilde{\phi}(u)\right)}\\\
\tilde{\phi}(u)=\phi(u)+{\bf\delta_{0}}.{\bf
u_{r}}\frac{d\phi_{0}}{du}-\psi_{1}(u)=\phi_{0}(u)+\epsilon\tilde{\psi}(u)\\\
\tilde{\psi}(u)=\psi(u)+{\bf\delta}.{\bf
u_{r}}\frac{d\phi_{0}}{du}-\psi_{1}(u)\\\
\nabla\psi_{1}={\bf\delta}\end{cases}$ (13)

According to Eq. (3) we have, $f_{n}=\frac{1}{n!}\frac{\partial\psi}{\partial
r}$. The notation $\tilde{\psi}$ represent the field $\psi$ once the effect of
the offset is taken into account. As a consequence for the offset fields
$f_{1}$ and $\frac{df_{0}}{d\theta}$ using Eq. (13) we obtain,

$\begin{cases}\tilde{f}_{1}=\left[\frac{\partial\tilde{\psi}}{\partial
u}\right]_{u=1}=f_{1}+{\bf\delta_{0}.u_{r}}\left(-1+\left[\frac{d^{2}\psi_{0}}{dr^{2}}\right]_{u=1}\right)=f_{1}-\kappa_{2}{\bf\delta_{0}.u_{r}}\\\
\frac{d\tilde{f}_{0}}{d\theta}=\left[\frac{\partial\tilde{\psi}}{\partial\theta}\right]_{u=1}=\frac{d\tilde{f}_{0}}{d\theta}+{\bf\delta_{0}}.\frac{\partial{\bf
u_{r}}}{\partial\theta}\left[\frac{d\phi_{0}}{dr}\right]_{r=1}-{\bf\delta_{0}.u_{\theta}}=\frac{df_{0}}{d\theta}\end{cases}$
(14)

Note that in Eq (14) we used $\nabla\phi_{1}={\bf\delta}$ to obtain the radial
and angular derivatives of $\phi_{1}$ (see Eq. (13)). We also used the fact
that $\left[\frac{d\phi_{0}}{dr}\right]_{r=1}$ which is the consequence of the
critical condition for the unperturbed potential at the scale radius
$R_{E}=1$. Here we used the variable $\kappa_{2}$ which was already used in
the initial presentation of the singular perturbative theory, see Alard
(2007). It is interesting to note that to the first order Eq. (14) indicates
that the field $\frac{df_{0}}{d\theta}$ is unaffected by the shift
$\delta_{0}$. The field $f_{1}$ itself is affected by the shift.

###  3.2 The choice of the referential.

In general the Fourier expansion of the perturbative fields at first order
contain four independent terms. Let’s define these first order terms as
${\bf\alpha}$ and ${\bf\beta}$ for $f_{1}$ and $\frac{df_{0}}{d\theta}$. We
have also to include the source impact parameter $\bf r_{0}$ (see Eq. 4). The
first order Fourier of a perturbative fields $F$ is represented by the suffix
$F_{1}$,

$\begin{cases}\left[\tilde{f}_{1}\right]_{1}={\bf(\alpha+r_{0})u_{r}}\\\
\left[\frac{d\tilde{f}_{0}}{d\theta}\right]_{1}={\bf(\beta+r_{0})u_{\theta}}\end{cases}$
(15)

We will consider that in our optimal referential these first order
(${\bf\alpha,\beta}$) Fourier terms are not zero but are the smallest
possible. If we consider a slight change in referential by introducing a small
shift $\bf\lambda$ by using Eq. 14 we obtain,

$\begin{cases}\left[\tilde{f}_{1}\right]_{1}={\bf(-\kappa_{2}\lambda+r_{0}+\alpha)u_{r}}\\\
\left[\frac{d\tilde{f}_{0}}{d\theta}\right]_{1}={\bf(\beta+r_{0})u_{\theta}}\end{cases}$
(16)

It is interesting to note that Eq. (16) implies a fundamental degeneracy in
the choice of the referential. In effect one could choose $\lambda$ such that,

$-\kappa_{2}\lambda+\alpha=\beta$ (17)

Since in the optimal referential model $|\bf\beta|$ is small then by direct
consequence of Eq. (17) the coefficients at first order for $f_{1}$ are small
two. Thus by this small shift we found another referential where the first
order terms are small. The basic idea behind this choice is that most dark
matter halos are not too distant from an isothermal profile. Furthermore it is
simple to show that for slightly non-isothermal halos the fields are
conjugates for first order Fourier terms.

####  3.2.1 Expression of the perturbative fields for slightly non-isothermal
potentials.

For a centered circular isothermal potential both perturbative fields are
null. When the same potential is slightly off-centered and shifted by
${\bf\lambda}$ at the first order the potential will be,

$\phi=|{\bf r-\lambda}|\simeq|{\bf r}|+{\bf\lambda u_{r}}$ (18)

The corresponding perturbative fields can be estimated by using
$f_{1}=\left[\frac{d\phi}{dr}\right]_{r=1}-1+{\bf r_{0}u_{r}}$ and
$\frac{df_{0}}{d\theta}=\left[\frac{d\phi}{d\theta}\right]_{r=1}+{\bf
r_{0}u_{\theta}}$. Here $r_{0}$ is the source impact parameter (see Eq. 4).

$\begin{cases}\left[\tilde{f}_{1}\right]_{1}={\bf(\lambda+r_{0})u_{r}}\\\
\left[\frac{d\tilde{f}_{0}}{d\theta}\right]_{1}={\bf(\lambda+r_{0})u_{\theta}}\end{cases}$
(19)

If we consider a set of slightly of centered isothermal with corresponding
weights potential Eq. (19) indicates that the contribution to the perturbative
fields will be the weighted sum of of the of the impact parameters of each
component. Thus for such slightly perturbed isothermal potentials we see that
the perturbative fields will have conjugate expressions and depends only on
two terms. As a consequence we will choose the conjugate field model in the
optimal referential. This is a reasonable choice for dark matter halos in this
degenerate problem. In any case as long as the shift from the optimal
referential is small and of the order of the perturbation size, the specific
choice of the referential should not matter and should not affect the final
result.

###  3.3 Estimating a first guess for the solution.

The process of building the initial guess is de-composed in two steps, first
an estimation of the Einstein ring size and center, and second an estimation
of the first guess of perturbative fields based on a round source model.

####  3.3.1 Estimating the position and size of the Einstein ring.

The first step is make some estimation of the radius and center of the
Einstein circle associated to the lens. The most direct method is to fit a
circle to the arc system. To fit a circle to the data we will consider a model
where the arc system is a slightly off-centered in the current coordinate
system. As a first guess for the center of the circle we may consider the
center of gravity the arc system or more directly the center of the image. The
actual choice of the guess is not very important as the search for the center
and radius of the circle converge quickly even for quite a distant initial
guess. For a circle of radius $R_{c}$ off-centered by $\bf r_{c}$ to the first
order in $\frac{r_{c}}{R_{c}}$ the circle equation reads,

$r\simeq R_{c}+{\bf
r_{c}u_{r}}=R_{c}+\left(x_{c}\cos\theta+y_{c}\sin\theta\right)$ (20)

The first step of our processing is to reduce the image of the arc to a set of
positions defined by the radius and corresponding angle. The processing starts
by choosing a threshold. In the calculations we will consider only the pixels
with values above this threshold. In practice the angular domain is divided
into a number of bins $N_{\theta}$ and for each bin the mean radius
$R_{\theta}$ is estimated. Then the First order Fourier series defined by Eq.
(20) is fitted to the set of positions $(R_{\theta},\theta)$. This fitting
procedure estimates the parameters $r_{0}$, and $\bf r_{c}$. Note that first
we estimate the parameter $\bf\alpha$ in Eq. (20) and then we infer $\bf
r_{c}$. Once the estimation of the parameters has been performed the
referential is shifted according to the estimation of $\bf r_{c}$. The the
process is iterated, a new estimation is performed and the center is refined.
This fitting process is very efficient and converges quickly. The result is a
first estimation of the Einstein radius and a first estimation of the lens
center. The circle fitting procedure converges when the circle is centered in
the referential system. When the method converges and the circle is centered
in the referential system the first order Fourier terms in Eq. (20) are zero.
For a circular source model the field $f_{1}$ represents the mean position of
the images (see Eq. 5 and Sec. 2.1). As a consequence for a simple round
source model the first order Fourier terms of the fields $f_{1}$ have to be
zero in the final referential after the circle fitting procedure has
converged. In a system shifted by a vector $\delta_{0}$ from the optimal
referential Eq. (14) and Eq. (19) implies that,

$\begin{cases}\left[\tilde{f}_{1}\right]_{1}={\left(\bf\lambda+r_{0}+\kappa_{2}\delta_{0}\right)u_{r}}={\left(\bf\bar{r}_{0}-\kappa_{2}\delta_{0}\right)u_{r}}\\\
\left[\frac{d\tilde{f}_{0}}{d\theta}\right]_{1}={\bf(\lambda+r_{0})u_{\theta}}={\bf\bar{r}_{0}u_{\theta}}\end{cases}$
(21)

In Eq. (21) the degenerate term ${\bf\lambda+r_{0}}$ was re-defined by
introducing the ${\bf\bar{r}_{0}=\lambda+r_{0}}$. Note that in the referential
defined by the circle fitting procedure the first order terms have to be zero,
which by using Eq. (21) is equivalent to,

$\bar{r}_{0}=\kappa_{2}\delta_{0}$ (22)

Our choice for the centering of the arcs is based on an isothermal halo
conjugated model (See Sec. 3.2.1).For an isothermal halo we have,
$\kappa_{2}=1$, $\bf\bar{r}_{0}=r_{0}$, Eq. (22) implies that
$\delta_{0}=r_{0}$ which means that the shift $\delta_{0}$ corresponds to the
source impact parameter. As a consequence for an isothermal potential the
circle centering procedure implies that the referential is centered on the
source. As illustrated by the rotation curves the dark matter halos are not
far from isothermal thus we the referential will be close to a source centered
referential. An interesting point is that at first order the field
$\frac{d\tilde{f}_{0}}{d\theta}$ is un-affected by the shift and consequently
allows a direct estimation of the pseudo impact parameter $\bf\bar{r}_{0}$
(see Eq. (22)).

####  3.3.2 Estimating a first guess for the perturbative fields.

A first rough reconstruction of the perturbative fields in this lens centered
initial referential will be performed. In the continuation the estimation of
the angular bin mean positions and width will be all renormalized by the
estimated Einstein radius. In the circular source model the field $f_{1}$ is
simply the image mean position while the field $\frac{df_{0}}{d\theta}$ is
related to the image width $\Delta$,

$\Delta=\sqrt{R_{0}^{2}-\frac{df_{0}}{d\theta}^{2}}$ (23)

The reconstruction of the field $f_{1}$ is obtained directly by fitting a
Fourier series to the mean position $R_{\theta}$ of the binned data at each
angular bin position where an image is available. For the field
$\frac{df_{0}}{d\theta}$ the data to fit will be the image width $\Delta$ (Eq.
23) for each angular bin. An estimation of $\Delta$ will be obtained by
counting the number of pixel in each angular bin above the threshold. The
number of pixels above the threshold is a direct estimation of the bin surface
from which the effective width $\Delta$ can be easily estimated. In practice
the number of pixel above a threshold is an integer quantity and as a
consequence may not be accurate for a small number of pixels. To improve the
accuracy of this estimate the number of pixels just below the threshold and
just above the threshold are estimated (see Eq. (26)). An interpolation of
these two nearby pixel counts provide an accurate estimation of the pixels
number at the threshold for each angular bin $C_{\theta}$.

$\Delta=\frac{N_{\theta}C_{\theta}}{2\pi R_{c}^{2}}$ (24)

It is in principle possible to evaluate the field $\frac{df_{0}}{d\theta}$ by
relating the numerical data of Eq. (24) with equation Eq. (23). However in
practice a direct evaluation of $\frac{df_{0}}{d\theta}$ using Eq’s (23) and
(24) is problematic for two reasons. The first reason is that due to the
square in Eq. (23) the sign of $\frac{df_{0}}{d\theta}$ is unknown. Each time
the function $\frac{df_{0}}{d\theta}$ cross the zero line the sign may or may
not change. The field $\frac{df_{0}}{d\theta}$ is zero near the center of each
image, thus if there are many images the uncertainty on the sign means that
the number of potential solution to explore is large. And furthermore the
second issue is that the management of the noise is very delicate due to the
square root in Eq. (23). An optimal estimation of $\frac{df_{0}}{d\theta}$
requires a fitting procedure based on a least-square model. The Fourier series
are the simplest and most direct non parametric reconstruction of the
perturbative fields since they are directly related to the multipole expansion
on the Einstein circle (see Eq. (6)). To construct the initial guess it is
important to reduce the complexity of the solution to a minimum and thus to
consider an expansion in Fourier series of the fields of lower order. A
typical number for the order of the Fourier $N_{F}$ expansion is $N_{F}=2$ or
$N_{F}=3$ depending on the complexity of the solution. It would be very
unusual and quite impractical to reach $N_{F}=4$ for the initial guess.
Another point is the existence of a possible mass-sheet degeneracy related to
the scaling coefficient $\kappa_{2}$ (see Eq. (3)). In order to cancel the
effect of the re-scaling coefficient $\kappa_{2}$ it is useful to introduce
the correlation coefficient which is scale free rather than a least-square
estimate. As a consequence the estimation of $\frac{df_{0}}{d\theta}$ will be
performed by maximizing the correlation coefficient between the two expression
of $\Delta$, the analytical expression in Eq. (23) and the numerical data of
Eq. (24). The order of the Fourier series is low and the cost of computing the
correlation coefficient for the array of binned data is low since in practice
the number of bins is less than one hundred. As a consequence it is definitely
practical to explore the parameter space to estimate a model for
$\frac{df_{0}}{d\theta}$. The field $\frac{df_{0}}{d\theta}$ has no zero order
term since it is a derivative, thus at order 2 the number of coefficients is 4
and 6 at order 3. The parameter space is explored by generating a series of
random numbers in a typical range. The range to investigate is by nature of a
fraction of unity. In practice the range at order 1 is about 0.5, and decrease
by a scale factor corresponding to the order at higher order. The assumption
of decreasing order is justified by the fact that the power of the Fourier
series at some higher order will be close to zero. Thus between the first
order and this higher order a decrease is observed and we make the choice of
the simplest model, a linear decrease. In practice this choice is notation
very important since we start with a very broad range at order 1 which is much
larger than practical cases. This model of linear decrease can thus be
considered as an upper bound. At order 2 the exploration of the parameter
space is performed using $10^{6}$ set of random number. The corresponding
computing time is only of about a few seconds on a basic PC for 75 angular
bins. Considering the low computing cost the number of random sets can be
easily increased if higher order Fourier solutions have to be explored.

####  3.3.3 Illustration using a numerical simulation.

The application of this method to a numerical simulation of a cusp
configuration for an elliptical NFW potential is presented in Fig. (1) (upper
left panel). The details of this numerical simulation are available in Sec.
(3.6.1). The different steps of the reconstruction are illustrated in the
upper right panel with the circle fitting procedure, and in the lower panel
with the first guess for the perturbative fields based on the round source
model.

Figure 1: Upper left panel: the simulated image corresponding to simulation
$A_{1}$. Upper right panel: the reconstructed image. Lower left panel: the
first guess estimation for the field $\frac{df_{0}}{d\theta}$ (black cross)
superimposed to the true solution(red cross). lower right panel: the first
guess estimation for the field $f_{1}$ (black cross) superimposed to the true
solution(red cross).

###  3.4 Refining the guess and converging to an optimal solution.

Once a a first approximate solution based on the circular source model has
been obtained for the fields the problem is to refine this solution and find
the optimal solution. In practice the approximate solution will be used as a
starting guess for a full minimization in the multi-parameter space. However
the first thing to consider in the search for an optimal solution is to define
an optimal referential. An optimal referential should minimize the non-
circular parts of the lens potential and as a consequence should be closely
centered on the lens. However the circle fitting procedure implies that the
first guess is defined in a source centered referential (see Sec. (3.3.1) and
Eq’s (21) and (22)). Consequently the first step should be to move from the
source centered referential where we defined the initial guess to a lens
centered referential. This task is facilitated by the fact that the field
$\frac{df_{0}}{d\theta}$ is unchanged by a small shift (see Eq. (14)). When
expressed in a referential shifted by $\bf\delta_{0}$ the first order term of
$\frac{df_{0}}{d\theta}$ are precisely the pseudo-impact parameter
$\bf{\bar{r}_{0}}$ (see Eq. (21)). Thus from the guess solution for
$\frac{df_{0}}{d\theta}$ we can estimate the pseudo-impact parameter
$\bar{\bf}r_{0}$. A problem is that the guess solution is estimated from the
circular solution of Eq. (5). In this equation the square of
$\frac{df_{0}}{d\theta}$ appear which mean that a degeneracy on sign of
$\frac{df_{0}}{d\theta}$ is present. This means that the sign of
$\bf\bar{r}_{0}$ is unknown as well as the higher order Fourier terms of the
field $\frac{df_{0}}{d\theta}$. In practice it mean that the two solutions
corresponding to the two signs of $\frac{df_{0}}{d\theta}$ have to be
explored. The circular source solution is insensitive to the sign of
$\frac{df_{0}}{d\theta}$ but this is not the case for a non-circular source.
As a consequence the $\chi_{2}$ will have to be evaluated for each sign of
$\frac{df_{0}}{d\theta}$ and the minimal $\chi_{2}$ will point point to the
proper sign and break the degeneracy.

####  3.4.1 Reconstructing the lens and source by minimizing the $\chi_{2}$.

####  3.4.2 The parameters of the lens reconstruction.

The parameters that will have to be refined are the following, the Einstein
radius $R_{C}$, the center of the unperturbed circular potential $\bf r_{c}$,
the pseudo impact parameter $\bar{\bf}r_{0}$, and the higher order Fourier
terms ($N>1$) of the perturbative fields. Note that the first order Fourier
terms of the perturbative fields are conjugated (see Eq. (21) here for
$\bf\delta_{0}=0$), and as a consequence there are only two independent terms.
The first order terms for $f_{1}$ are directly inferred from the first order
terms of $\frac{df_{0}}{d\theta}$ by conjugation.

####  3.4.3 The source reconstruction procedure.

The source is represented by the linear combination of 2D Gaussian polynomials
(see Alard & Lupton (1998)). The 2D amplitude in the source plane represented
by the coordinates $(x_{s},y_{S})$ is,

$I_{S}=\sum_{i=0,N}\sum_{j=0,N-i}a_{i,j}\exp(-\left(\frac{x_{S}^{2}+y_{S}^{2}}{r_{0}^{2}}\right)x^{i}y^{j}$
(25)

The modelisation of complex sources which are frequently encountered in
gravitational lenses requires a higher order for the expansion in Eq. (25). In
practice to avoid the source-lens degeneracy the initial source expansion
should be of low order, but should be increased to reach a final higher order
with typically $N=8$ or $N=10$. An initial estimation of $r_{0}$ is obtained
by estimating the residual to the fit of a circle to the arc (see Sec 3.3.1).
Once a first fit of the lens has been performed $r_{0}$ is evaluated directly
by computing the size of the actual source solution. This process of refining
$r_{0}$ is iterated for each steps of the fitting procedure.

####  3.4.4 Fitting the image and estimating the $\chi^{2}$.

For a given model of the perturbative fields to Fourier order the use of the
lens equation allows (Eq. (3)) allows to transport each component of the
source model (Eq. 25) to the lens plane. Each corresponding image of the
source component is then convolved with a Gaussian PSF model. The set of
convolved images is then fitted to the arcs system images using the linear
least-square method. This fitting procedure allows the evaluation of the
$\chi_{2}$ corresponding to the specific model for the perturbative fields.

###  3.5 Iterating the procedure.

Starting from the initial circular source guess we will make a search for the
minimal $\chi_{2}$ solution. The minimization of the $\chi_{2}$ in the
parameter space of the lens (see Sec. (3.4.2) ) is performed using the simplex
method (Nelder & Mead (2006)).For each set of paramaters the $\chi_{2}$ is
evaluated according to the procedure described in Sec. (3.4.4). Due to the
degeneracy on the sign of $\frac{df_{0}}{d\theta}$ the search for the minimal
$\chi_{2}$ is conducted for both negative and positive signs. The choice
between the 2 solutions with different signs is made by selecting the solution
with minimal $\chi_{2}$.

####  3.5.1 Avoiding the source lens degeneracy problem.

The source lens degeneracy may be present due to the interaction and
degeneracy between higher order source and lens terms. As shown in Alard
(2017) a solution to this issue is to consider the minimal solutions. These
minimal solutions are among the class minimal $\chi_{2}$ solutions those with
the lowest possible order in the Fourier expansion of the fields. It is
demonstrated in Alard (2017) that these minimal solutions can be constructed
by starting the search for the solution to the lowest possible order. A
typical number to reconstruct the guess solution is $N=2$ or $N=3$ for the
Fourier expansion of the fields. This will be also the order at which the
solution is reconstructed for the first step of the refinement. Once the sign
degeneracy is broken the order of the solution is increased progressively to
reach the highest order for the solution (typically $N=4$ or $N=5$). Note that
in similarity to the reconstruction process of the guess solution the size of
the simplex corresponding to the higher order terms is of decreasing
amplitude. Here the same rule applies, with a decrease in the amplitude
proportional to order. It is possible to relax this constraint for some
specific case, however it is better to try using this modulation of the higher
order terms first.

###  3.6 Application of the method to a set of numerical simulations.

In this paragraph we will demonstrate the refinement of the guess solution
using a set of 4 numerical simulations.

####  3.6.1 Description of the numerical simulation.

The four numerical simulations corresponds to a cusp like configuration for 4
different potential types. The different potential models are a second order
lens potential (elliptical) a fourth order lens potential (quadratic) and an
external shear. The combination of the elliptical and quadratic potential with
the external shear give the four potential models (see Table 1) and the four
corresponding arcs simulations $A_{i=1,n}$.

Arcs Simulation | Potential Fourier order | External Shear
---|---|---
$A_{1}$ | 2 | no
$A_{2}$ | 2 | yes
$A_{3}$ | 4 | no
$A_{4}$ | 4 | yes
Table 1: The four types of simulations.

####  3.6.2 The potential models in the simulations

The potential models are based on the NFW potential described in Eq. (28) in
(Alard (2007)), and (Bartelmann (1996)). The iso-contour of this potential are
slightly distorted from circularity. The distortion is expressed using Fourier
series to order $N$. The coefficients of the Fourier series representing the
distortion are represented by $\eta_{i}$. The lowest order corresponds to an
elliptical distortion, $N=2$ (2 coefficients), and the highest order in this
simulation corresponds to $N=4$ (quadratic).

$\begin{cases}\phi\propto\frac{1}{2}\log^{2}\left(\frac{x}{2}\right)-2\operatorname{arctanh}^{2}\left(\sqrt{\frac{1-x}{1+x}}\right)\\\
x=r\left(1+\sum_{k}\eta_{[2(k-2)]}\cos(k\theta)+\eta_{[2(k-2)+1]}\sin(k\theta)\right)\end{cases}$
(26)

The shear is external and corresponds to the effect of mass situated outside
the Einstein circle. The shear is related to order 2 terms in the potential
and as a consequence to order 2 Fourier series terms in the potential
expansion. Using Eq. (8) for an external shear of order 2 he generic
expression for the fields is,

$\begin{cases}\left[f_{1}\right]_{shear}=S_{0}\cos 2\theta-S_{1}\sin
2\theta\\\ \left[\frac{df_{0}}{d\theta}\right]_{shear}=-S_{1}\cos 2\theta-
S_{0}\sin 2\theta\end{cases}$ (27)

For the two numerical simulation with shear we have $S_{0}=S_{1}=0.05$.

####  3.6.3 Building the simulations.

For each of the potential models $A_{i}$ the image of the source is
reconstructed using ray-tracing. The source in this simulation is complex and
quite distant from a round source (see Fig. (2) lower panel for an image of
the source contours). The image of the source was convolved with a Gaussian
PSF with a FWMH of 2.5 pixels. The final image was produced by adding Poisson
noise to the convolved image.

####  3.6.4 Elliptical second order potential with no shear.

This simulation corresponds to case $A_{1}$ in Table (1). The process of
building this guess solution was described in Sec. (3.3). The resulting guess
solution for this simulation is illustrated in Fig. (1). Starting from this
guess solution the solution is refined in a series of steps with each time an
increase of the order for the Fourier series expansion of the potential. The
refinement operates according to the procedure described in Sec. (3.4.1). The
refinement process starts by exploring the two solutions corresponding to each
sign of the field $\frac{df_{0}}{d\theta}$ (see Sec (3.5)). In this first
iteration the Fourier order of the potential is 2. From the two solutions the
solution with minimal $\chi_{2}$ is retained. This solution is refined by
increasing the Fourier order of the potential to 3 and subsequently to 4 to
avoid the source degeneracy problem (see Sec. (3.5.1)). The refinement process
is stopped at order $N=4$. Not that for minimality reason the reconstruction
could have been stopped at order $N=2$ in this case. However in order to
provide a comparison of the different lens reconstructions in a self
consistent scheme all lens are reconstructed to order 4. The reconstruction of
the lens to a higher is also interesting since it demonstrates the stability
of the solution at higher order. All solutions were built using an eight order
Gaussian-polynomial solution for the source (See Sec. (LABEL:Sec_source)).The
distortion vector corresponding to simulation $A_{1}$ is,

$\eta=[-0.05,0,0,0,0,0]$ (28)

The final result is presented in Fig.’s (2) and (3).

Figure 2: Upper left panel: a comparison between the reconstructed fields
(black line) and the true solution (red line) Upper right panel: The amplitude
of the Fourier series expansion of the fields as a function of the order. In
both figures the field $\frac{df_{0}}{d\theta}$ is represented as a dotted
line, while the $f_{1}$ is represented by a continuous line. The lower left
panel represents the true source while the lower right panels represents the
reconstructed source. Figure 3: Upper left panel: the image of the arc system
corresponding to simulation $A_{1}$. Upper right panel: the reconstructed
image obtained after the final iteration. Lower left panel: subtraction
between the true and reconstructed image. Lower right panel: an histogram the
noise normalized by the theoretical noise expectation. Note that the histogram
(red line) is quite close to the Gaussian expectation (black line).

####  3.6.5 Elliptical second order potential with shear.

The simulation $A_{2}$ is similar to simulation $A_{1}$ except for the
presence of an external shear (see Eq. (27)). The construction of the guess
solution and of the final refined solution follows the same process that was
already described for simulation $A_{1}$. The final result of the refinement
process is presented in Fig.’s (4) and (5).

Figure 4: Same as a Fig. (2) but this time for simulation $A_{2}$. Figure 5:
Same as a Fig. (3) but this time for simulation $A_{2}$.

####  3.6.6 Fourth order potential with no shear.

Here simulation $A_{3}$ corresponds to an NFW potential with a fourth order
Fourier distortion of the iso-contours . The corresponding distortion vector
$\eta_{i}$ i define in Eq. (26). Here again the construction of the guess
solution and of the final refined solution follows the same process that was
already described for simulation $A_{1}$. For simulation $A_{3}$ we have,

$\eta=[-0.05,0,0.01,0,0,-0.005]$ (29)

The final result of the refinement process is presented in Fig.’s (6) and (7).

Figure 6: Same as a Fig. (2) but this time for simulation $A_{3}$. Figure 7:
Same as a Fig. (3) but this time for simulation $A_{2}$.

####  3.6.7 Fourth order potential with shear.

Simulation $A_{4}$ corresponds to simulation $A_{3}$ with the addition of an
external shear(see Eq. (27)). Here again the construction of the guess
solution and of the final refined solution follows the same process that was
already described for simulation $A_{1}$. The final result of the refinement
process is presented in Fig.’s (8) and (9).

Figure 8: Same as a Fig. (2) but this time for simulation $A_{4}$. Figure 9:
Same as a Fig. (3) but this time for simulation $A_{4}$.

####  3.6.8 Result of the simulations.

It is clear form the results of the simulations see Fig’s (3), (5), (7) and
(9) that all configurations and even the more complex were properly
reconstructed using the automated method presented in this paper. No
assumptions were made about the arcs themselves and the procedure produced the
lens and source modelling directly from the images. All results were produced
using the same universal singular perturbative approach which suggests that an
automated analysis of a large data set will now become possible with this
method. Interestingly the universality of the modelling will also render
possible a statistical analysis of large data sets.

##  4 Conclusion.

In this paper systems of gravitational arcs of increasing complexity were
analysed and successfully reconstructed. Even in the case of a fourth order
lens potential and a strong external shear an accurate reconstruction of the
potential and source was achieved. The same complex highly non symmetrical
source model was used for all the simulations. This analysis demonstrates that
the method developed in this paper has the ability to perform automated
reconstructions of gravitational arcs systems. The reconstructions are direct
due to the nature of the singular perturbative formalism and as a consequence
insensitive to degeneracy issues. Considering the rapidly increasing number of
good quality images of gravitational arcs systems this package of programs is
of particular interest. The method offers the the possibility to reconstruct
quickly a large number of arcs using a common universal model. The
universality of the potential modelling allows the reconstruction of
statistical quantities for a large number of arcs. This possibility is of
particular interest when the data from incoming high quality sky surveys like
Euclide will be available.

## Data Availability

No data sets were generated or analyzed during the current study.

## Acknowledgements

## References

  * Agnello etal. (2013) Agnello, A., Auger, M. W., Evans, N. W., 2013, MNRAS, 429L, 35
  * Alard & Lupton (1998) Alard, C., Lupton, R., 1998, ApJ, 503,325
  * Alard (2007) Alard, C., 2007, MNRAS, 382, L58
  * Alard (2008) Alard C., 2008, MNRAS, 388, 375
  * Alard (2009) Alard C., 2009, A&A, 506, 609
  * Alard (2010) Alard C., 2010, A&A, 513, 39
  * Alard (2017) Alard C., 2017, MNRAS, 467, 1997
  * Alard (2017) Alard C., 2023, MNRAS, 521, 4896
  * Bartelmann (1996) Bartelmann, M., 1996, A&A,313,697
  * Brewer etal. (2016) Brewer, Brendon J.; Huijser, David; Lewis, Geraint F., 2016, MNRAS, 455, 1819
  * Meneghetti etal. (2006) Meneghetti, M., Bartelmann, M., Moscardini, L., 2003, MNRAS, 340, 105
  * Nelder & Mead (2006) Nelder, J., & Mead, R. 1965, Comp. J., 7, 308
  * Saha & Williams (2006) Saha, P., Williams, L., 2006, ApJ, 653, 936
  * Warren & Dye (2003) Warren, S., & Dye, S. 2003, ApJ, 590, 673