#  Faster Convergence for Transformer Fine-tuning with Line Search Methods
††thanks: We gratefully acknowledge the funding by the German Federal Ministry
of Economic Affairs and Energy (01MK20007E). The work was partially conducted
within the KI-Starter project “Robustness and Generalization in Training Deep
Neural Networks” funded by the Ministry of Culture and Science Nordrhein-
Westfalen, Germany

1st Philip Kenneweg  AG Machine Learning
University Bielefeld
Bielefeld, Germany
0000-0002-7097-173X 2nd Leonardo Galli  MIP Chair
RWTH Aachen University
Aachen, Germany
0000-0002-8045-7101 3rd Tristan Kenneweg  AG Machine Learning
University Bielefeld
Bielefeld, Germany
0000-0001-8213-9396 4th Barbara Hammer  AG Machine Learning
University Bielefeld
Bielefeld, Germany
0000-0002-0935-5591

###### Abstract

Recent works have shown that line search methods greatly increase performance
of traditional stochastic gradient descent methods on a variety of datasets
and architectures [1, 2]. In this work we succeed in extending line search
methods to the novel and highly popular Transformer architecture and dataset
domains in natural language processing. More specifically, we combine the
Armijo line search with the Adam optimizer and extend it by subdividing the
networks architecture into sensible units and perform the line search
separately on these local units. Our optimization method outperforms the
traditional Adam optimizer and achieves significant performance improvements
for small data sets or small training budgets, while performing equal or
better for other tested cases. Our work is publicly available as a python
package, which provides a hyperparameter-free pytorch optimizer that is
compatible with arbitrary network architectures.

###### Index Terms:

Transformer, Stochastic Line Search, Optimization, BERT, NLP

##  I INTRODUCTION

In modern machine learning, there is a variety of optimization algorithms.
Figuring out which one is the best and matching the learning rate or learning
rate schedule to the specific problem can require a lot of expert knowledge
and computing power. In particular, the current state of the art is to treat
the learning rate as a hyperparameter and train the network until the value
that yields the best performance is found, i.e., hyperparameter tuning. To
simplify and significantly accelerate this process, a recent branch of the
deep learning research [3, 4, 5, 1] has suggested to reintroduce line search
methods for selecting the step size. These methods are widely adopted in the
optimization community since their introduction in [6]. In particular, they
automatically find an adaptive learning rate by calculating the loss at
different points along the (gradient) direction and selecting an appropriate
step size. Using line search can save time and resources by eliminating the
need for manual optimization of the learning rate, thus eliminating costly
hyperparameter tuning.

As traditional line search requires multiple forward passes per gradient
update, a more efficient approach is needed. In [1], a Stochastic Line Search
(SLS) has been combined with a smart re-initialization of the step size to
alleviate the need for multiple forward passes for every step. This approach
was shown to outperform a variety of other optimization methods, such as
Stochastic Gradient Descent (SGD), SGD with Nesterov Momentum and many others
on tasks such as matrix factorization as well as image classification for
small networks and datasets.

It is however unclear if stochastic line search methods work well for larger
networks and different architectures such as Transformers [7] as well as in
challenging domains such as Natural Language Processing (NLP). In this work we
adapt and expand upon this active topic of research in the NLP domain.

In addition, given the deep structure of recent networks, it is natural to ask
if it would be beneficial to train each layer at a different speed, i.e.
learning rate, instead of using a single value for the whole network. For this
reason, we propose to modify SLS [1] and apply it in a layer-wise manner in
order to specialize the step size for separate network components. This novel
idea aims at localizing the line search in order to capture the loss variation
due to the update of a single network component.

The Adam [8] optimizer has been shown to outperform SGD in the training of
Transformers [9] by a large margin. In this work, we propose to use SLS [2] in
combination with Adam to train Transformers or any other architecture. In
particular, we will compare the following optimization schemas on the problem
of fine-tuning Transformers:

  * •

Adam optimizer with warm starting and cosine decay (ADAM).

  * •

Armijo line search combined with stochastic gradient descent (SGDSLS).

  * •

Armijo line search combined with the Adam optimizer (ADAMSLS).

  * •

Per-layer-ADAMSLS, as a novel optimization method (PLASLS).

We test all optimizers fine-tuning BERT[10] on the Glue[11] dataset
collection, which is widely used to evaluate common natural language
processing skills. We find that SGDSLS does not work well in this scenario. In
contrast, ADAMSLS and PLASLS perform better than ADAM, especially if they are
trained on smaller data sets, as is often the case for retraining language
models [10]. To make our work easy to reproduce and easy to use, we implement
all methods as pytorch optimizers. ADAMSLS can be used as an optimizer, that
does not require to manually set the learning rate. For PLASLS the network
parameters need to be manually split into components which are optimized
separately.

In the next Section, we will review the related literature, in Section III we
will describe our method. In Section IV we will describe our experimental
approach and in Section V we will show and discuss our results. In Section VI
we will go into details on some choices for PLASLS. Finally Section VII
contains our conclusion.

##  II RELATED WORK

The optimization of deep neural networks has been a central topic of research
in the field of machine learning. A variety of techniques and optimizers have
been proposed, including but not limited to SGD [12], Adagrad[13], RMSprop[14]
and Adam[8]. One common challenge in optimizing neural networks is the
selection of appropriate step sizes, which can have a significant impact on
the convergence rate and final performance of a model [15].

Line search methods are a popular approach for selecting step sizes in
optimization algorithms. These methods involve iteratively adjusting the step
size based on the curvature of the loss function, in order to ensure that the
objective is decreasing at each iteration. However, traditional line search
methods can be computationally expensive, particularly for large neural
networks with many layers.

In this work we particularly build upon [1]. In this paper, the authors show
the theoretical proofs of convergence for the stochastic Armijo line search
method in the case of strongly convex, convex and non convex functions.
Furthermore, they provide solutions to some practical problems of line search
methods. In particular, they no longer start the search from a fixed and
possibly high initial step size, but instead from the last selected step size.
To allow the step to also increase, they double it every $m/b$ steps, where
$b$ is the minibatch size and $m$ is a constant that does normally not need to
be tuned, see Eq. 2. Additionally, they show empirical results on the image
datasets MNIST[16], CIFAR10 [17] and CIFAR100[17] using the ResNet[18] and
DenseNet[19] architecture, as well as on a variety of convex problems.

Important advantages of the Armijo line search compared to other optimization
methods are:

  * •

no hyperparameter tuning of the learning rate

  * •

faster convergence rates

  * •

better generalization

However, their work does not show whether Armijo line search can outperform
optimizers such as Adam on more complex tasks and architectures. To
investigate this, we focus on natural language processing with Transformers,
as this is a critical area of recent development and high complexity where
other promising optimization methods have been unable to improve upon the
baseline set by the Adam optimizer. Additionally, we investigate whether
splitting the networks global step size into locally optimized step sizes is a
viable option to further improve the performance of line search methods.

Recent work has shown that Transformers are very sensitive to the learning
rate and learning rate schedule schedule they are trained with[20, 21]. To
remedy this, various approaches like RADAM[21] or warm starting have been
proposed. We show that our approach is able to train these highly sensitive
architectures well.

Other recent work includes optimizing the learning rate simultaneously with
the network weights [22]. The authors showed that this yields good results on
classical image datasets, but does not seem to speed up convergence.

Further related work includes [23], which focuses on layer specific learning
rates, [24] which studies the scaling limit of SGD in the high dimensional
regime and [9] which studies why Adam is so effective at training the
Transformer architecture.

Overall, the optimization of neural networks and transformer-based models
remains an active area of research, and the development of efficient line
search methods is an important contribution to this field.

##  III METHODS

The classical stochastic Armijo line search [1] is designed to set a maximum
step size for all network parameters $w_{k}$ at iteration $k$ (Eq. 1). In this
section, together with the classical SGD direction, we will include the more
complex Adam [8] direction. Moreover, we instead suggest to select a different
step size for each different ”component” of the network and evaluate the loss
decrease w.r.t. it.

We define the following notation: The loss function is denoted by $f(w)$.
$||\cdot||$ denotes the Euclidean norm and $\nabla f$ denotes the gradient of
$f$. Given the iteration counter $k$, $f_{k}$ and $\nabla f_{k}$ denote the
mini-batch function and its mini-batch gradient.

###  III-A Armijo Line Search

The Armijo line search criterion is defined in [1] as:

$f_{k}(w_{k}+\eta_{k}d_{k})\leq f_{k}(w_{k})-c\cdot\eta_{k}||\nabla
f_{k}(w_{k})||^{2},$ (1)

where $d_{k}$ is the direction (e.g., $d_{k}=-\nabla f_{k}(w_{k})$ in case of
SGD), $c>0$ is a constant (commonly fixed to be $0.1$ [1, 9, 2]). Condition 1
is practically obtained by employing a backtracking procedure, i.e., starting
with an initial step-size $\eta^{0}_{k}$ and iteratively decreasing it by a
constant factor $\delta\in(0,1)$ until Eq. 1 is satisfied. As described in
Section II, the following tweak was introduced in [1] to reduce the amount of
backtracks, but avoid a monotonically decreasing step size

$\eta_{k}=\eta_{k-1}\cdot 2^{b/m}.$ (2)

###  III-B Including Adams Update Step in SLS (ADAMSLS)

In case of SGD, the direction $d_{k}$ is the negative mini-batch stochastic
gradient.

$d_{k}=-\nabla f_{k}(w_{k})$

Adam’s direction [8] can be defined as

$\begin{split}g_{k+1}&=\nabla f_{ik}(w_{k})\\\ m_{k+1}&=\beta_{1}\cdot
m_{k}+(1-\beta_{1})\cdot g_{k+1}\\\ v_{k+1}&=\beta_{2}\cdot
v_{k}+(1-\beta_{2})\cdot g_{k+1}^{2}\\\
\hat{m}_{k+1}&=m_{k}/(1-\beta^{k}_{1})\\\
\hat{v}_{k+1}&=v_{k}/(1-\beta^{k}_{2})\\\
d_{k}&=-\hat{m}_{k+1}/(\sqrt{\hat{v}_{k+1}}+\epsilon)\end{split}$ (3)

Adam combines a momentum-based approach together with a step size correction
built upon the gradients variance. In the training of Transformers, these
modifications have been shown to be important enhancements over the simpler
SGD[24]. The weight update rule is generally defined as

$w_{k+1}=w_{k}+\eta_{k}d_{k}.$ (4)

To apply the Armijo line search on the Adam optimizer we now use the direction
$d_{k}$ defined in Eq. 3, but with momentum $\beta_{1}=0$. Additionally, the
gradient norm term $||\nabla f_{k}(w_{k})||^{2}$ is replaced by the scaled
gradient norm $\frac{||\nabla
f_{k}(w_{k})||^{2}}{\sqrt{\hat{v}_{k+1}}+\epsilon}$ resulting in Eq. 5.

$f_{k}(w_{k}+\eta_{k}d_{k})\leq f_{k}(w_{k})-c\cdot\eta_{k}\frac{||\nabla
f_{k}(w_{k})||^{2}}{\sqrt{\hat{v}_{k+1}}+\epsilon},$ (5)

###  III-C Layer Wise Line Search (PLASLS)

The key idea of layer wise line search is to split up the network parameters
into $L$ arbitrary sub-units, such as provided by network layers. In
particular, we split $w_{k}$, $d_{k}$ and the gradient $\nabla f_{k}(w_{k})$
into their components and rewrite $\eta_{k}$ as a vector:

$\begin{split}w_{k}&=(w_{k}^{(1)},\dots,w_{k}^{(l)},\dots,w_{k}^{(L)})\\\
\nabla f_{k}(w_{k})&=\left({\nabla f_{k}(w_{k})}^{(1)},\dots,{\nabla
f_{k}(w_{k})}^{(l)},\dots,{\nabla f_{k}(w_{k})}^{(L)}\right)\\\
d_{k}&=(d_{k}^{(1)},\dots,d_{k}^{(l)},\dots,d_{k}^{(L)})\\\
\eta_{k}&=(\eta_{k}^{(1)},\dots,\eta_{k}^{(l)},\dots,\eta_{k}^{(L)}).\end{split}$

At this point, we adapt the Armijo line search to involve only elements of the
sub-unit $(l)$. However, the loss $f$ is not partially separable [25], i.e.,
it can only be computed on the whole set of parameters $w_{k}$. Therefore, we
define gradient components $\bar{d}_{k,l}$ and rewrite Eq. 1 as follows

$\begin{split}f_{k}(w_{k}+\eta_{k}^{(l)}\bar{d}_{k,l})&\leq
f_{k}(w_{k})-c\cdot\eta_{k}^{(l)}||{\nabla f_{k}(w_{k})}^{(l)}||^{2},\\\
\bar{d}_{k,l}&:=(0,\dots,d_{k}^{(l)},\dots,0)\end{split}$ (6)

This not only reduces the computational costs of the line search, but also
localizes the method to capture the variation of the loss only w.r.t. the
$l$-th sub-unit. At each iteration $k$, we only apply a line search procedure
on a single sub-unit to yield $\eta_{k}^{(l)}$, while all the other step sizes
remains unchanged. In other words, if $l$ is the sub-unit selected at
iteration $k$ the update can be rewritten as

$\begin{split}w_{k+1}&=w_{k}+\eta_{k,l}\odot d_{k},\\\
\eta_{k,l}&:=(\eta_{k-1}^{(1)},\dots,\eta_{k}^{(l)},\dots,\eta_{k-1}^{(L)}),\end{split}$
(7)

where $\odot$ is the Hadamard product, or component-wise product. Thus the
resulting computational complexity is very similar to the original
implementation, while allowing different step sizes for different network
components.

#### Splitting

In the case of the Transformer[7] architecture or more specifically the
BERT[10] architecture which is encoder only, we split the network into 10
different components. BERT has 12 encoder layers we divide the different
layers into 8 different parts ($n_{k}^{(2)}$ affects layer 1, $n_{k}^{(3)}$
affects layers 2,3, $n_{k}^{(4)}$ affect layer 4, $n_{k}^{(5)}$ affect layers
5,6 …). $n_{k}^{(1)}$ affects the embedding weights of the network.
$n_{k}^{(10)}$ affects the single dense layer which is appended for a specific
task. This is the optimal split found as shown by the experiments in Section
VI.

#### First Step

Since the step size for each sub-unit is not changed during each iteration, a
large initial step size $\eta_{0}$ could lead to divergence problems. To make
our method completely automatic, we first perform a single step of SLS on the
whole network to select a common first initial step size $\eta_{0}$. In
particular, we choose a big initial step size $\eta_{0}=0.1$ and let the first
line search find a good replacement for it. In fact, the step size yielded by
the line search is a upper bound for the inverse of the Lipschitz constant of
$\nabla f(w)$ (see Lemma 1 in [1]) and this provides us with a good
approximation of the complexity of the problem at hand (at least in the
starting point). From this value, we then start subdividing the step sizes per
sub-unit.

#### Step Size Merging

If we look closely at the resulting step sizes yielded by our methods, we
observe that they might converge to very low values (e.g., Figure 1). This
also happens without the layer wise splitting, but less frequently. As a
result, the affected part of the network is no longer able to learn. This can
be desirable in some cases, but seems problematic overall.

Figure 1: Exemplary problematic step size of component 9 of the network during
a single training run on the QNLI dataset. The step size starts out in the
order of $10^{-4}$ but is lowered to values below $10^{-50}$.

In the specific case of layer-wise line search optimizer, a first solution for
this problem would be to merge the network parts based on their current step
size.

Algorithm 1 algorithm for merging network components

step sizes $n_{k}^{(l)}$, threshold $\lambda$

find the smallest step size $n_{k}^{(s)}\in n_{k}^{(l)}$

if $n_{k}^{(s)}\leq\lambda$ then

find the second smallest step size $n_{k}^{(s2)}\in n_{k}^{(l)}$

merge network components $s$ and $s2$

$n_{k}^{(new)}=(n_{k}^{(s)}+n_{k}^{(s2)})/2$

end if

In Algorithm 1, we detect if a network components step size $n_{k}^{(l)}$ is
below a certain threshold $\lambda$. If this is the case, this network
component is merged with the network component with the second lowest step
size. The new step size for this combined network component is the average of
both step sizes.

Figure 2: Different merging thresholds on one epoch of the MNLI dataset.
Standard error is indicated around each line.

In Figure 2, the results of different merging thresholds $\lambda$ are
visualized, they seem to indicate that no merging performs worse than even the
smaller merging thresholds of for example $\lambda=10^{-12}$. Consequently,
all experiments in Section IV for PLASLS are conducted using the merging
threshold $\lambda=10^{-12}$. With this threshold only very little automatic
merging occurs. On the smaller datasets (all $small$ datasets, MRPC) automatic
merging does not occur, while on the larger datasets (SST2, QNLI, MNLI) on
average about 3 of the 10 components are merged during a full 5 epoch training
run. We strongly suspect that numerical cancellation errors are the cause of
very low step sizes, we therefore plan to investigate this topic in more
details in the future.

##  IV EXPERIMENTAL APPROACH

In this section, we detail our experimental design to investigate the
performance of the proposed optimizers. We utilize the Huggingface library
[26] for implementation and the pre-trained Bert model ’bert-base-uncased’ for
all experiments.

###  IV-A Candidates

As a baseline comparison we evaluate the Adam optimizer with a peak learning
rate of 2e-5 and a cosine decay with warm starting for 10% of the total
training time, henceforth referenced to as $ADAM$. These values are taken from
the original paper [10] and present good values for a variety of
classification tasks, including the Glue[11] tasks upon which we are
evaluating.

As another baseline, we use Armijo line search as described in [1] with
stochastic gradient descent, we further refer to this approach as $SGDSLS$.

Next, we use the Armijo line search described in [1] and combine it with the
Adam optimizer[8] as described in Section III to obtain the approach further
referenced to as $ADAMSLS$.

The last option we are testing is the layer wise line search as described in
Section III combined with the Adam optimizer further referenced to as
$PLASLS$.

###  IV-B Implementation Details

The metaparameter choices for all experiments are as follows:

  * •

All models are trained for 5 epochs.

  * •

The pooling operation used in all experiments is [CLS].

  * •

Batch size used for training is 32.

  * •

The Adam optimizer with betas (0.9,0.999) and epsilon $10^{-8}$ is used.

  * •

The maximum sequence length is set to 256 tokens.

  * •

All models are trained 5 times and their mean metrics and standard error are
reported.

###  IV-C Datasets

We consider a common scenario in natural language processing, where a large
pre-trained language model is fine-tuned on a small dataset. The Glue dataset
by Wang et al. [11] is a collection of various popular classification tasks in
NLP, and it is widely used to evaluate common natural language processing
capabilites. All datasets used are the version provided by tensorflow-datasets
4.0.1.

More specifically, we use the datasets Stanford Sentiment Treebank _SST2_ ,
Microsoft Research Paraphrase Corpus _MRPC_ , Stanford Question Answering
Dataset _QNLI_ , and Multi-Genre Natural Language Inference Corpus _MNLI_.

We evaluate all approaches on the $small$ and full size datasets. The $small$
datasets are the same as previously described, except that the size of the
training dataset has been reduced to 500 randomly drawn samples. This scaling
enables us to judge the capability of the optimizer in different dataset size
regimes, as especially in real world scenarios small dataset sizes are a
common occurrence.

(a) MNLI small

(b) MRPC small

(c) QNLI small

(d) SST2 small

Figure 3: The loss curves of the experiments on the small dataset with
standard error indicated around each line. In all experiments SGDSLS performs
the worst, followed by ADAM. PLASLS and ADAMSLS do not perform very different.
In the SST2 and QNLI experiment PLASLS performs best, while in the MNLI
experiment ADAMSLS performs best. In the MRPC experiments both perform about
the same.

##  V EXPERIMENTAL RESULTS

In this section, we will describe the results of our experiments on the Glue
dataset collection. We compare the 4 candidates as described in Section III

  * •

ADAM

  * •

SGDSLS

  * •

ADAMSLS

  * •

PLASLS

All accuracies displayed are the accuracies on the validation sets. The losses
displayed are the losses calculated on the training sets, smoothed with
exponential moving average. The colored areas around each line indicate the
standard error of each experiment. We display the accuracies and losses during
the training period in Figures 3,4,5,6. The Tables I and II are displaying the
accuracies averaged over all datasets used, as is commonly done for the Glue
dataset collection.

(a) MNLI small

(b) MRPC small

(c) QNLI small

(d) SST2 small

Figure 4: The accuracy curves of the experiments on the small dataset with
standard error indicated around each line, starting after the first epoch. In
all experiments SGDSLS performs the worst, followed by ADAM. PLASLS and
ADAMSLS do not perform very different. In the MNLI and MRPC experiment ADAMSLS
performs best. In the SST2 and QNLI experiments ADAMSLS and PLASLS perform
about the same.

###  V-A Small Experiments

In the small size experiments (see Figure 3 and Figure 4) we observe that the
final loss of ADAMSLS and PLASLS is at least one order of magnitude lower
compared to ADAM and SGDSLS. This is also reflected in the accuracy metric
where PLASLS and ADAMSLS perform significantly better than ADAM or SGDSLS.
Small losses on the training data and high accuracies on the validation data
are highly correlated and no overfitting is observed.

In Table I, we see the average performance over all datasets and runs at the
end of training. ADAMSLS and PLASLS clearly outperform ADAM and SGDSLS by
about 3%. Note that this accuracy is taken after 5 epochs, earlier in training
even larger advantages for ADAMSLS and PLASLS can be observed.

TABLE I: Average classification accuracies, for the $small$ datasets. Best performing method is marked in bold.  method | ADAM | SGDSLS | ADAMSLS | PLASLS
---|---|---|---|---
$accuracy$ | 0.6875 | 0.6927 | 0.7250 | 0.7165

###  V-B Full Size Experiments

(a) MNLI

(b) MRPC

(c) QNLI

(d) SST2

Figure 5: The loss curves of the experiments on the full size dataset with
standard error indicated around each line. Overall SGDSLS clearly performs
worst. In the SST2 experiments PLASLS fails to converge to a very low loss. In
the MRPC experiment we can see that ADAMSLS and PLASLS perform better
initialy, but ADAM performs about the same in the end.

(a) MNLI

(b) MRPC

(c) QNLI

(d) SST2

Figure 6: The accuracy curves of the experiments on the full size dataset with standard error indicated around each line, starting after the first epoch. Accuracy was calculated on the validation data. Overall we can observe that SGDSLS performs the worst. ADAMSLS is overall the best performing optimizer.  TABLE II: Average classification accuracies, for the fullsize datasets. Best performing method is marked in bold.  method | ADAM | SGDSLS | ADAMSLS | PLASLS
---|---|---|---|---
$accuracy$ | 0.8745 | 0.8714 | 0.8830 | 0.8779

In the full size experiments, we observe that the loss starts to decrease
quicker with ADAMSLS and PLASLS. After about 10k steps ADAM outperforms the
other options on the loss metric. On the accuracy metric we can see that
ADAMSLS and PLASLS perform better in the beginning, but only slightly better
at the end of training.

In Table II we see the final average performance over all datasets and runs.
ADAMSLS and PLASLS very slightly outperform ADAM and SGDSLS by about 1%. It
should be noted that these accuracies are measured after 5 epochs. Earlier in
training, greater benefits can be observed for ADAMSLS and PLASLS.

###  V-C Discussion

Overall we see that with smaller training sets or shorter training runs,
ADAMSLS and PLASLS perform significantly better than ADAM or SGDSLS. This
performance difference vanishes for longer training runs, here ADAM seems to
perform better on the loss metric but still performs slightly worse in
accuracy.

Interestingly, we can see in Figure 7 that the average step size of a well
hyperparameter tuned ADAM and the automatically calculated ADAMSLS and PLASLS
are similar during long periods of training. This demonstrates, the ability of
ADAMSLS and PLASLS to result in very similar step size regimes to the tuned
BERT steps size for the Adam optimizer even though they are not given an
initial step size and thusly do not need any hyperparameter tuning.

Figure 7: Average step size of all layers during training runs on the QNLI
dataset.

##  VI ABLATION STUDIES

In this section we will take a short look at which parts of our PLASLS
approach have the most impact on performance, as well as different options for
multiple elements of the described algorithm which we tested, but did not go
into detail in the other sections.

###  VI-A Batch Size

During our experiments we noticed that the batch size had a large influence on
the stability of the training. In this section we compare different batch
sizes.

Figure 8: Loss on the sst2small dataset for different batch size options for
PLASLS. Standard error is visualized around each line. The X-Axis denotes the
number of samples processed.

In Figure 8, we see that PLASLS performance improves for larger batch sizes
during training. Furthermore, the loss curves get a lot smoother with higher
batch sizes. During our experiments we choose a batch size of 32 as it was the
highest supported batch size on our available hardware.

###  VI-B Split Options

Several possibilities to split a Transformer based network exist. In this
section we look at a split by:

  * •

layer

  * •

query, key and value

We tried splitting the network in 1,4,7,10 components as well as by query,key
and value (QKV). If we split the network in one component, ADAMSLS is
equivalent to PLASLS. The loss curves for different layer configurations can
be seen in Figure 9. More layers result in gradually better performance for
the SST2small dataset. Splitting by QKV does not seem to improve performance.

Figure 9: Loss on the SST2small dataset for different split options for PLASLS
and ADAMSLS as reference. Standard error is visualized around each line.

##  VII CONCLUSION

The type of optimizer, learning rate and learning rate schedule is a critical
choice for every machine learning pipeline. We are the first in the literature
to evaluate the combination of Adam and line search methods on transformer
fine tuning(ADAMSLS) and present a new variant of the line search optimizer
PLASLS. Both optimizers we presented perform equal or better than state-of-the
art baselines, especially on small data sets or short training runs in the
domain of natural language processing for the ubiquitous Transformer
architecture. We recommend using our ADAMSLS PyTorch implementation as best
practice for the task of Transformer fine tuning.

For further work it would be interesting to combine the faster initial
convergence rate of PLASLS or ADAMSLS with the long-term convergence of ADAM.
Either by finding a more stochasticly stable replacement for the Armijo
criterion, or by manually switching to the Adam optimizer after a certain
amount of steps. By combining/improving the optimizer thusly, we could get an
optimizer which would be even more widely applicable.

The source code is open-source and free (MIT licensed) software and available
at
https://github.com/TheMody/Faster-Convergence-for-Transformer-Fine-tuning-
with-Line-Search-Methods.

## References

  * [1] S. Vaswani, A. Mishkin, I. Laradji, M. Schmidt, G. Gidel, and S. Lacoste-Julien, “Painless stochastic gradient: Interpolation, line-search, and convergence rates,” NIPS’19: Proceedings of the 33rd International Conference on Neural Information Processing Systems, 2019.
  * [2] S. Vaswani, F. Kunstner, I. H. Laradji, S. Y. Meng, M. Schmidt, and S. Lacoste-Julien, “Adaptive gradient methods converge faster with over-parameterization (and you can do a line-search),” CoRR, vol. abs/2006.06835, 2020.
  * [3] M. Mahsereci and P. Hennig, “Probabilistic line searches for stochastic optimization,” Advances in neural information processing systems, vol. 28, 2015.
  * [4] R. Bollapragada, J. Nocedal, D. Mudigere, H.-J. Shi, and P. T. P. Tang, “A progressive batching l-bfgs method for machine learning,” in International Conference on Machine Learning, pp. 620–629, PMLR, 2018.
  * [5] C. Paquette and K. Scheinberg, “A stochastic line search method with expected complexity analysis,” SIAM Journal on Optimization, vol. 30, no. 1, pp. 349–376, 2020.
  * [6] L. Armijo, “Minimization of functions having lipschitz continuous first partial derivatives,” Pacific Journal of mathematics, vol. 16, no. 1, pp. 1–3, 1966.
  * [7] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. u. Kaiser, and I. Polosukhin, “Attention is all you need,” in Advances in Neural Information Processing Systems (I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, eds.), vol. 30, Curran Associates, Inc., 2017.
  * [8] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings (Y. Bengio and Y. LeCun, eds.), 2015.
  * [9] F. Kunstner, J. Chen, J. W. Lavington, and M. Schmidt, “Noise is not the main factor behind the gap between sgd and adam on transformers, but sign descent might be,” in The Eleventh International Conference on Learning Representations, 2023.
  * [10] J. Devlin, M. Chang, K. Lee, and K. Toutanova, “BERT: pre-training of deep bidirectional transformers for language understanding,” in Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, NAACL-HLT 2019, Minneapolis, MN, USA, June 2-7, 2019, Volume 1 (Long and Short Papers) (J. Burstein, C. Doran, and T. Solorio, eds.), pp. 4171–4186, Association for Computational Linguistics, 2019.
  * [11] A. Wang, A. Singh, J. Michael, F. Hill, O. Levy, and S. Bowman, “GLUE: A multi-task benchmark and analysis platform for natural language understanding,” in Proceedings of the 2018 EMNLP Workshop BlackboxNLP: Analyzing and Interpreting Neural Networks for NLP, (Brussels, Belgium), pp. 353–355, Association for Computational Linguistics, Nov. 2018.
  * [12] H. Robbins and S. Monro, “A stochastic approximation method,” The annals of mathematical statistics, pp. 400–407, 1951.
  * [13] J. Duchi, E. Hazan, and Y. Singer, “Adaptive subgradient methods for online learning and stochastic optimization,” Journal of Machine Learning Research, vol. 12, no. 61, pp. 2121–2159, 2011.
  * [14] G. H. with Nitish Srivastava Kevin Swersky, “Lecture notes neural networks for machine learning,” 2014.
  * [15] K. Nar and S. S. Sastry, “Step size matters in deep learning,” in Proceedings of the 32nd International Conference on Neural Information Processing Systems, NIPS’18, (Red Hook, NY, USA), p. 3440–3448, Curran Associates Inc., 2018.
  * [16] Y. LeCun and C. Cortes, “MNIST handwritten digit database,” 2010.
  * [17] A. Krizhevsky, “Learning multiple layers of features from tiny images,” pp. 32–33, 2009.
  * [18] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 770–778, 2016.
  * [19] G. Huang, Z. Liu, L. Van Der Maaten, and K. Q. Weinberger, “Densely connected convolutional networks,” in 2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 2261–2269, 2017.
  * [20] P. Kenneweg, A. Schulz, S. Schröder, and B. Hammer, “Intelligent learning rate distribution to reduce catastrophic forgetting in transformers,” in Intelligent Data Engineering and Automated Learning – IDEAL 2022 (H. Yin, D. Camacho, and P. Tino, eds.), (Cham), pp. 252–261, Springer International Publishing, 2022.
  * [21] L. Liu, H. Jiang, P. He, W. Chen, X. Liu, J. Gao, and J. Han, “On the variance of the adaptive learning rate and beyond,” in International Conference on Learning Representations, 2020.
  * [22] K. Chandra, A. Xie, J. Ragan-Kelley, and E. Meijer, “Gradient descent: The ultimate optimizer,” in Advances in Neural Information Processing Systems (A. H. Oh, A. Agarwal, D. Belgrave, and K. Cho, eds.), 2022.
  * [23] B. Singh, S. De, Y. Zhang, T. Goldstein, and G. Taylor, “Layer-specific adaptive learning rates for deep networks,” 10 2015.
  * [24] G. B. Arous, R. Gheissari, and A. Jagannath, “High-dimensional limit theorems for SGD: Effective dynamics and critical scaling,” in Advances in Neural Information Processing Systems (A. H. Oh, A. Agarwal, D. Belgrave, and K. Cho, eds.), 2022.
  * [25] L. Galli, A. Galligari, and M. Sciandrone, “A unified convergence framework for nonmonotone inexact decomposition methods,” Computational Optimization and Applications, vol. 75, no. 1, pp. 113–144, 2020.
  * [26] T. Wolf, L. Debut, V. Sanh, J. Chaumond, C. Delangue, A. Moi, P. Cistac, T. Rault, R. Louf, M. Funtowicz, J. Davison, S. Shleifer, P. von Platen, C. Ma, Y. Jernite, J. Plu, C. Xu, T. L. Scao, S. Gugger, M. Drame, Q. Lhoest, and A. M. Rush, “Huggingface’s transformers: State-of-the-art natural language processing,” 2019.


# Modulation of the solar microwave emission by sausage oscillations

Elena G. Kupriyanova,1 Tatyana I. Kaltman,2 and Alexey A. Kuznetsov3
1Central Astronomical Observatory at Pulkovo of RAS, Saint Petersburg, 196140,
Russia
2Saint Petersburg Branch, Special Astrophysical Observatory of RAS, Saint
Petersburg, 196140, Russia
3Institute of Solar-Terrestrial Physics SB RAS, Irkutsk, 664033, Russia
E-mail<EMAIL_ADDRESS>(EGK)

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The modulation of the microwave emission intensity from a flaring loop by a
standing linear sausage fast magnetoacoustic wave is considered in terms of a
straight plasma slab with the perpendicular Epstein profile of the plasma
density, penetrated by a magnetic field. The emission is of the
gyrosynchrotron (GS) nature, and is caused by mildly relativistic electrons
which occupy a layer in the oscillating slab, i.e., the emitting and
oscillating volumes do not coincide. It is shown that the microwave response
to the linear sausage wave is highly non-linear. The degree of the non-
linearity, defined as a ratio of the Fourier power of the second harmonic to
the Fourier power of the principal harmonic, is found to depend on the
combination of the width of the GS source and the viewing angle, and is
different in the optically thick and optically thin parts of the microwave
spectrum. This effect could be considered as a potential tool for diagnostics
of the transverse scales of the regions filled in by the accelerated
electrons.

###### keywords:

Sun: radio radiation – (magnetohydrodynamics) MHD – radiation mechanisms: non-
thermal – Sun: activity

††pubyear: 2022††pagerange: Modulation of the solar microwave emission by
sausage oscillations–Modulation of the solar microwave emission by sausage
oscillations

##  1 Introduction

Solar flares are one of the most powerful phenomena in the solar atmosphere.
Physical processes operating in flares at kinetic and magnetohydrodynamic
(MHD) scales, such as the charged particle acceleration and magnetic
reconnection, are fundamental processes of plasma astrophysics and remain
subject to intensive studies. An important method for probing the plasma in
flaring sites is MHD seismology, based upon oscillatory processes detected in
flares (e.g. Kupriyanova et al., 2020; Zimovets et al., 2021). Among the
oscillations which are most commonly observed in flaring regions are standing
fast magnetoacoustic waves of the sausage symmetry, i.e. the axisymmetric
oscillatory motions across the magnetic field, characterised by variations of
the plasma density and absolute value of the magnetic field (see, e.g., Li et
al., 2020, for recent comprehensive reviews). As the typical periods of
sausage oscillations are from a few seconds to several tens of seconds (e.g.,
Nakariakov & Melnikov, 2009), their detection requires the use of instruments
with high time resolution. In particular, sausage oscillations are often
observed as the modulation of coherent and non-coherent radio and microwave
emissions. For example, sausage oscillations have been used for the
interpretation of the time variations of the gyrosynchrotron (GS) emission
intensity from spatially resolved coronal loops Nakariakov et al. (2003);
Melnikov et al. (2005), and simultaneous oscillations of GS and soft X-ray
emission along a loop (Inglis et al., 2008), the wiggling of zebra-pattern
lanes (Yu et al., 2013), and the precipitation rate of the nonthermal
electrons at the opposite footpoints of a loop (Nakariakov et al., 2018).

Sausage oscillations are used as seismological probes of perpendicular
profiles of the density and the absolute value and twisting of the magnetic
field in flaring plasma structures, such as plasma loops (e.g., Li et al.,
2020). Both regular and random perpendicular profiles could be addressed (e.g.
Pascoe et al., 2007; Hornsey et al., 2014; Chen et al., 2015; Lopin & Nagorny,
2015b). Typically, the analysed plasma structures are modelled as a plasma
slab or a cylinder, surrounded by a plasma with different properties (see,
e.g. Nakariakov et al., 2012; Hornsey et al., 2014), respectively. The
consideration of sausage oscillations in a loop with a varying cross-section
demonstrated that the slab or cylinder models are adequate. For example, for
the fundamental parallel harmonic, the cross-section radius increase near the
loop top by a factor of 2 causes the decrease in the oscillation period of
about 5% only (Pascoe et al., 2009).

Despite intensive theoretical and observational studies of sausage
oscillations, the link of the theory and observational outcomes remains non-
trivial. One of the difficulties is the effect of the line-of-sight (LoS)
integration in the optically thin regime. Let us illustrate it for the case of
the EUV or soft X-ray emissions. On one hand, a sausage oscillation is
characterised by the perturbations of the plasma density. One the other hand,
if the LoS is perpendicular to the oscillating plasma structure, the plasma is
displaced along the LoS. As the observed emission intensity is proportional to
an integral of the squared density along the LoS, the observed intensity does
not change much in different phases of the oscillation. This effect was first
demonstrated by Gruszecki et al. (2012), and then confirmed by comprehensive
forward modelling (Antolin & Van Doorsselaere, 2013). A similar problem
appears in the radio observations (see, e.g. Mossessian & Fleishman, 2012),
and has been addressed by forward modelling of microwave observables by
Reznikova et al. (2015); Kuznetsov et al. (2015).

In contrast with the modelling of the thermal emission which is determined by
the variation of the plasma density along the LoS, various kinds of non-
thermal radio emissions are also determined by the magnetic field non-thermal
electrons, LoS angle to the field, and their distributions along the LoS. In
the majority of forward modelling studies, the spatial structure of non-
thermal electrons has been considered to coincide with the oscillatory
perturbations of the plasma density and magnetic field by a sausage wave
Reznikova et al. (2015); Kuznetsov et al. (2015). However, this assumption
could be incorrect, as the volumes occupied by non-thermal emission and
involved in the sausage oscillation are determined by different physical
processes. The perpendicular structure of the sausage wave perturbation is
determined by parameters of wave-guiding plasma non-uniformity (e.g., Li et
al., 2020). The magnetic flux tube filled in with non-thermal electrons is
determined by the acceleration mechanism. For example, a coronal structure,
which is observed in microwaves as a single loop, appears as a filamented
structure in EUV (e.g., Zimovets et al., 2013). This structure can oscillate
as a whole, being perturbed by a sausage wave, while only one its filaments
may contain energetic electrons emitting the microwaves.

The mismatch of the oscillating and emitting volumes could provide us with
important seismological constrains on the enigmatic non-thermal electron
acceleration processes. The microwave emission could be considered as a tool
for diagnostics of the low-amplitude perturbation of the plasma magnetic by
the sausage wave because of the high sensitivity and the non-linear dependence
of the intensity of the microwave emission on the both magnitude and direction
of the magnetic field and on the parameters of the emitting electrons (e.g.,
Dulk & Marsh, 1982), resulting in the non-linear microwave response to the
linear sausage wave. The aim of this paper is to study how the degree of this
non-linearity depends on both the width of the emitting volume and the LoS
angle to the magnetic field.

The paper is organised as follows. The model and basic equations are described
in Section 2. The source of the microwave emission and details of the
calculation of the microwave emission are presented in Section 3, and details
of obtaining the microwave light curves are presented in Section 4. Section 5
presents the principal results of the simulations, which are summarized and
discussed in Section 6.

##  2 Plasma slab and MHD wave

We consider a transversely inhomogeneous (along the $x$-axis) plasma slab
stretched along a uniform magnetic field $B_{0}$ directed along $z$-axis. The
slab is perturbed by a symmetric (sausage) fast magnetoacoustic MHD mode of
the slab. The system is uniform in the $y$ direction. If we place an
infinitely remote observer in the $xz$ plane, then we may ignore the $y$
direction and simulate the MHD wave in the two-dimensional slab. So, the
components of the unperturbed magnetic field are
$B_{z0}=B_{0}=\textrm{const}$, $B_{x0}=0$.

The transverse profile of the plasma non-uniformity, shown by black curve in
panel (a) in Figure 1, is defined by the Epstein function (Cooper et al.,
2003)

$\rho_{0}(x)=\rho_{\mathrm{max}}\rm{sech}^{2}\it\frac{x}{w}+\rho_{\infty}$ (1)

with the characteristic width $w$ referred to as the slab width, and the
density contrast

$d=\frac{\rho_{\mathrm{\max}}+\rho_{\infty}}{\rho_{\infty}}.$ (2)

The system is symmetric relatively to $z$ axis (blue line in Figure 1) where
the plasma density is $\rho_{\mathrm{\max}}+\rho_{\infty}=\rho_{0}(x=0)$.
Plasma density at the point infinitely remote from the $z$ axis is
$\rho_{\infty}=\rho_{0}(x\to\infty)$.

Under conditions typical of the solar corona, where the plasma parameter
$\beta\ll 1$, the ideal MHD equations (e.g., Nakariakov et al., 1997) can be
linearised, and the perturbed plasma density and the perturbed magnetic field
components can be expressed via the transverse component of the perturbed
speed $\tilde{V}_{x}=\tilde{V}_{x}(x,z,t)$:

$\tilde{\rho}=-\int\frac{\partial(\rho_{0}\tilde{V_{x}})}{\partial x}dt,$ (3)
$\tilde{B_{x}}=B_{0}\int\frac{\partial\tilde{V_{x}}}{\partial
z}dt,\>\>\>\>\>\>\>\>\>\>\>\>\tilde{B_{z}}=-B_{0}\int\frac{\partial\tilde{V_{x}}}{\partial
x}dt.$ (4)

Since we consider the standing sausage wave, which is periodic in the $z$
direction with the longitudinal wave number $k$ and wavelength
$\lambda=2\pi/k$, the general expression of the transverse component
$\tilde{V}_{x}$ (see, e.g., Nakariakov & Roberts, 1995) can be simplified:

$\tilde{V}_{x}(x,z,t)=AU(x)\sin(kz)\cos(kV_{\mathrm{ph}}t).$ (5)

Here $A$ is the amplitude of the MHD perturbation normalized by the Alfvén
speed $C_{\mathrm{A0}}$ at the slab axis ($x=0$), $A\ll 1$, and
$V_{\mathrm{ph}}=\omega/k$ is the phase speed. Analytical expression for the
transverse structure $U(x)$ of the sausage wave

$U(x)=\frac{\sinh(x/w)}{\cosh^{\zeta}(x/w)},\>\>\>\>\>\>\>\>\>\zeta=\frac{|k|w}{C_{\mathrm{A\infty}}}\sqrt{C_{\mathrm{A\infty}}^{2}-V_{\mathrm{ph}}^{2}}+1,$

for the plasma density inhomogeneity defined by Equation (1) was obtained by
Cooper et al. (2003). The phase speed is defined by the dispersion equation

$\frac{|k|w}{C_{\mathrm{A0}}^{2}}(V_{\mathrm{ph}}^{2}-C_{\mathrm{A0}}^{2})-\frac{2}{|k|w}=\frac{3}{C_{\mathrm{A\infty}}}\sqrt{C_{\mathrm{A\infty}}^{2}-V_{\mathrm{ph}}^{2}},$

here $C_{\mathrm{A\infty}}$ is the Alfvén speed at the infinitely remote point
($x\to\infty$).

Following Equations (3)–(5), the perturbed values for plasma density and the
components of the magnetic field become

$\rho=\rho_{0}+\tilde{\rho},\>\>\>\>\>\>\>\>\>B_{z}=B_{0}+\tilde{B}_{z},\>\>\>\>\>\>\>\>\>B_{x}=\tilde{B}_{x}.$
(6)

Distributions of these parameters in the $xz$ plane at the moment of the
maximal perturbation are shown in panels (b), (c), (d) in Figure 1 where axes
$x$ and $z$ are normalised by the characteristic width of the plasma non-
uniformity $w$. In this paper, for simplicity, we consider strictly
monochromatic (i.e. strictly repetitive) decayless MHD wave; therefore, it is
sufficient to consider one oscillation period. The wavelength of the MHD wave
was chosen to be $\lambda=10w$.

Figure 1:  Panel (a): normalised transverse profile of the non-perturbed
thermal plasma density. The vertical short line indicates the characteristic
slab width $w$. Two rectangles illustrate homogeneous distributions of the
density of the accelerated electrons within both the narrow GS source
($|l/w|=0.5$, grey area) and the wide GS source ($|l/w|=3$, cross-hatched
area). Panel (b): perturbed plasma density $\rho$. Panel (c): perturbed
$B_{z}$ component of the magnetic field. Panel (d): perturbed $B_{x}$
component of the magnetic field. The darkest colour corresponds to the maximal
value, while the lightest colour — to the minimal value. The green point
$z_{1}$ in panel (b) denotes one of the sausage wave nodes and the blue point
$z_{2}$ — one of the anti-nodes. Different pairs of the same-coloured wavy
lines delimit the GS sources having different widths.

##  3 Microwave emission source

Thermal plasma in the slab is the source of thermal free-free radio emission.
In addition, if there is a population of electrons accelerated up to mildly
relativistic energies, such a system is the source of gyrosynchrotron
emission, which is usually associated with solar radio emission in the
microwave band. The GS source is populated with energetic electrons
isotropically distributed over pitch-angles and power-law distributed over
energy; the concentration of the electrons is assumed to be constant in the
unperturbed source which, according to the above-described approach, occupies
a part of the sausage-oscillating slab. The GS source is enclosed between
pairs of the magnetic field lines symmetric relatively to the $z$ axis. In the
unperturbed slab, the magnetic field lines are straight lines at the distances
of $x/w=l/w$ from the $z$ axis; thus $|l/w|$ is the half-width of the GS
source. We consider the GS sources both narrower ($|l/w|\leq 1$) and wider
($|l/w|>1$) than the characteristic slab width. Examples of the unperturbed
electron density distributions are shown in Figure 1 (panel (a)) by both grey
rectangle for the narrow GS source ($|l/w|=0.5$) and the cross-hatched
rectangle for the wide GS source ($|l/w|=3$). Hereinafter, we will omit the
modulus sign from the width of the GS source, for the simplicity.

The sausage wave perturbs the parameters of the slab and, therefore, modifies
the shape of the magnetic field lines and, thereby, curves the boundaries of
the GS source along the slab. Examples of the pairs of the symmetric perturbed
magnetic field lines, which are the outer boundaries of the GS sources, are
shown in Figure 1 (panels (b), (c), (d)) by the same-coloured wavy curves. The
local density of the non-thermal particles was assumed to be inversely
proportional to the local width of the GS source, to provide that the flux of
the accelerated electrons through the cross-section of the slab is conserved.
So, the modulation of the local distance between two symmetric magnetic field
lines leads to the modulation of the electron density along $z$ axis.
Modulation of the energy and pitch-angle distributions by the MHD oscillations
was neglected.

We use the Fast Gyrosynchrotron Codes (Fleishman & Kuznetsov, 2010; Kuznetsov
& Fleishman, 2021) to calculate the radio emission from the slab. As an input
data, these codes use the specified distributions of the emission source
parameters (such as the thermal plasma density and temperature, concentration
and distribution characteristics of the nonthermal electrons, and the magnetic
field strength and direction) along a line-of-sight. The codes include the
contributions of both the GS (nonthermal) and free-free (thermal) emission
mechanisms. However, we note that the thermal emission is several orders of
magnitude weaker than the GS emission, for the model parameters listed below.
The codes calculate the resulting radio emission by integrating the radiation
transfer equations for both the left-hand ($L$) and right-hand ($R$) polarised
components, which correspond to either ordinary or extraordinary
electromagnetic modes depending on the local magnetic field direction.

We consider the following parameters of the unperturbed slab: magnetic field
strength $B_{0}=200$ G, plasma temperature $T_{0}=10^{7}$ K, plasma density at
the $z$ axis $\rho_{\mathrm{max}}=5\times 10^{9}$ cm-3, plasma density outside
the slab $\rho_{\infty}=\rho_{\mathrm{max}}/(d-1)$, where the density contrast
$d=10$ corresponds to the plasma in flaring loops, the characteristic slab
width $w=3\times 10^{8}$ cm. The nonthermal electrons in the unperturbed GS
source are assumed to have the constant concentration of
$n_{\mathrm{b}}=10^{5}$ $\textrm{cm}^{-3}$, power-law energy distribution in
the range from 0.1 to 10 MeV with the spectral index of $\delta=3$, and
isotropic pitch-angle distribution. These parameters were chosen to provide
the GS emission spectral peak at around 3–5 GHz, which is a typical value in
solar flares. We consider the $z$ axis inclined relatively to the line-of-
sight by the viewing angle $\theta$ from 40∘ to 89∘.

##  4 Light curves

Intensity of the GS emission strongly and non-linearly depends on both the
local magnitude and direction of the magnetic field (see, for example,
Equation (1) in (Fleishman & Kuznetsov, 2010) or the empirical approximation
in (Dulk & Marsh, 1982)). Standing sausage wave modulates parameters of the GS
source which results in oscillations of the observed radio flux.

Application of the Fast GS Codes allows one to obtain the distribution of the
intensity ($I=R+L$) of the radio emission along the plane-of-sky which is
inclined relatively to the $z$ axis by the angle $90^{\circ}-\theta$. Running
the procedure for each time during the period of the sausage wave, one obtains
the temporal evolution of the intensity (light curve). We split the period of
the wave into 20 time segments, so, by definition, the period of the sausage
wave $P_{\mathrm{saus}}=20$ time units (t. u.). We consider the light curves
at two positions: $z_{1}$ which is a node of the sausage wave and $z_{2}$
which is an anti-node (see panel (b) in Figure 1). A few representative
examples of the light curves, smoothed with the cubic spline interpolation,
are presented in Figures 2–4 (upper row in each figure). The light curves are
coloured in green for $z_{1}$ and blue for $z_{2}$. We consider both the
optically thin emission at 17 GHz (Figures 2–3) and the optically thick
emission at 1.5 GHz (Figure 4). In each of these figures, the upper left panel
corresponds to a narrower GS source ($l/w=0.1$), while the upper right panel —
to a wider GS source ($l/w=3$ or $l/w=1.5$). Figures 2 and Figure 4 correspond
to $\theta=50^{\circ}$ while Figure 3 — to $\theta=70^{\circ}$.

##  5 Results

Figure 2: The upper row: the microwave light curves in the optically thin
regime ($f=17$ GHz) from the sausage-oscillating slab, both for the narrower
GS source, $l/w=0.1$ (panel (a)) and for the wider GS source, $l/w=3$ (panel
(b)). In each panel, the green and blue curves are obtained along line-of-
sight passing through the node $z_{1}$ and the anti-node $z_{2}$, respectively
(see panel (b) in Figure 1) at the viewing angle $\theta=50^{\circ}$. The
lower row: the normalized Fourier powers (asterisks) of the periods $P_{1}$
and $P_{2}$ found in the light curve associated with the anti-node (see blue
curves in the upper panels) for both the narrower GS source (panel (c)) and
the wider GS source (panel (d)). Figure 3: Same as in Figure 2, for the
optically thin regime ($f=17$ GHz) and viewing angle of $\theta=70^{\circ}$.

It is clearly seen in Figures 2–4 that the light curves deviate from a
sinusoidal shape, i.e., the response of the radio emission to the linear MHD
perturbation is nonlinear. Moreover, this deviation depends on the width of
the GS source $l/w$ and the viewing angle $\theta$. For example, for
$\theta=50^{\circ}$, the non-linearity is more pronounced for the narrower GS
source, which is manifested as a flattening of minima or maxima of the light
curves (e.g., blue curve in panel (a) of Figure 2). In contrast, for
$\theta=70^{\circ}$, the flattening and hence non-linearity is more pronounced
for the wider GS source (see both blue and green curves in panel (b) of Figure
3).

Figure 4: Same as in Figures 2–3, for the optically thick regime ($f=1.5$ GHz)
and viewing angle of $\theta=50^{\circ}$.

Figure 5: Distributions of the ratios ($\log_{10}\xi$) of the Fourier power of
the second harmonic with the period $P_{2}=10$ t. u. to the Fourier power of
the principal harmonic with the period $P_{1}=20$ t. u. relatively to the
width of the GS source ($l/w$) and viewing angles $\theta$ in both the
optically thin regime at 17 GHz (panel (a)) and the optically thick regime at
1.5 GHz (panel (b)). The lighter colours correspond to the higher values of
$\log_{10}\xi$. Lower row: the ratios $\xi$ for the selected viewing angles
$\theta$, marked by vertical lines in the upper panels, for optically thin
regime (panel (c)) and for optically thick regime (panel (d)).

To estimate the degree of the non-linearity, we calculated the modulation
depth

$\Delta=\frac{I_{max}-I_{min}}{I_{mean}},$

where $I_{max}$ and $I_{min}$ are the maximal and minimal intensities during
the MHD wave period, and $I_{mean}$ is the respective mean value. For the
narrower GS source, the modulation depths are of $\Delta\sim 100$–200% and
$\Delta\sim 60$–130% in the optically thin and optically thick regimes,
respectively. The high values of $\Delta$ are comparable with those usually
associated with quasi-periodic ejection of electrons (e.g., Kupriyanova et
al., 2016). The deep modulation can be explained by the location of the narrow
GS source near the $z$ axis where the oscillations magnitude of the waveguide
parameters is maximal.

For the wider GS source, the modulation depth is lower: $\Delta\sim 10$–20%
for both the optically thin and optically thick emissions. On the other hand,
the non-linear effect is much stronger: an extra peak or dip in the
lightcurves appears either at $t=P_{\mathrm{saus}}/4=5$ (as in the green
curves in panel (b) of Figure 2 and panel (a) of Figure 4 or as in the blue
curve in panel (b) of Figure 4) or at $t=3P_{\mathrm{saus}}/4=15$ (as in the
blue curves in panel (b) of Figure 2 and panel (a) of Figure 4). The amplitude
of the extra peak/dip varies from around 10% to 100% of the above-mentioned
modulation depth $\Delta$ for different combinations of $l/w$ and $\theta$.

Lower panels in Figures 2–4 demonstrate the ideal Fourier periodograms, where
one of the peaks corresponds to the principal period of the sausage wave at
$P_{1}=P_{\mathrm{saus}}=20$ t. u. The nonlinear response of the radio
emission mechanism to the MHD oscillation results in appearance of the second
harmonic with the period $P_{2}=10$ t. u. The principal period dominates in
most of the periodograms, while the second harmonic is less intensive.
However, for some combinations of $l/w$ and $\theta$, the Fourier power of the
second harmonic becomes higher; an example is presented in Figure 4 (panel
(d)).

We should stress that the second harmonic obtained in our study, which the
period exactly equals to the half-period of the principal harmonic,
$P_{2}=0.5P_{1}$, differs from the second harmonic of the eigen sausage mode,
which the period is greater than $0.5P_{1}$ because of the wave dispersion
(e.g., Nakariakov & Verwichte, 2005).

We characterize of the degree of non-linearity as a ratio ($\xi$) of the
Fourier power of the second harmonic to that of the principal harmonic.
Examples of the $\xi$ maps, as functions of the source width and the viewing
angle, are shown in Figure 5 for both the optically thin (panel (a)) and
optically thick (panel (b)) regimes. The maps correspond to the anti-node,
$z_{2}$ (see Figure 1); note that we plot $\log_{10}\xi$ values here. In the
optically thin regime, the $\xi$ values increase with the increasing of the
width $l/w$ for most angles $\theta$. Besides, there is a bright peculiarity
at the viewing angles around $\theta=63^{\circ}$, where the second harmonic
dominates over the principal harmonic. In the optically thick regime, the map
is more complicated.

Different characteristic properties of the ratio $\xi$ in both cases are also
clearly seen along the cross-sections of the maps at different selected angles
$\theta$ ($\theta=50$, 60, 70, and 80∘), marked by the vertical lines in two
upper panels of Figure 5. The cross-sections are presented in panels (c) and
(d) of Figure 5 for the optically thin and optically thick regimes,
respectively.

##  6 Summary and discussion

We performed simulations of the microwave radio response to a standing sausage
wave perturbing a zero-$\beta$ plasma slab, and, thereby, the source of the GS
emission, modulated by the sausage wave. The axis of the slab is directed
along the magnetic field, and the perpendicular profile of the plasma density
profile is smooth. The principal novel and important element which differs
this study from the previous ones is that the GS-emitting electrons fill only
a layer stretched along the axis of the slab, and localised in the
perpendicular direction, i.e., the oscillating and emitting volume do not
coincide. Our results confirm that the microwave response to a linear MHD wave
can be highly non-linear. This result was expected from, for example,
simplified expressions for the emission and absorption coefficients (Dulk &
Marsh, 1982) and was found previously in simulated microwave light curves
(Mossessian & Fleishman, 2012) in a simplified model. The principally new
result of our study is that the degree of the non-linearity depends on both
the relative width of the GS source inside the oscillating volume, and the
viewing angle. The non-linearity appears in the deviation of the microwave
intensity light curve from a sinusoidal shape of the sausage oscillation. This
results in appearing the second and higher harmonics of the monochromatic
sausage oscillation in Fourier spectra of the time-varying GS signal. Analysis
of the degree of non-linearity defined here as a ratio of the Fourier power of
the second harmonic to that of the principal harmonic allowed us to define
conditions where the non-linear effect is the strongest.

The effect found can be considered as a potential tool for the diagnostics of
the transverse scales of the regions containing accelerated electrons within
sausage-oscillating plasma structures. If we neglect the effect of non-thermal
electron diffusion across the magnetic field and take into account the
divergence of magnetic tubes in the cusp, this transverse size can be
considered as a characteristic perpendicular size of the acceleration region
and, therefore, as an upper limit for the perpendicular size of the magnetic
reconnection region. Note that previously this size was defined from numerical
simulations and was prescribed by certain MHD models. The results presented in
our paper open a possibility to estimate the spatial scales of the
reconnecting process independently, by analysing the QPPs in the microwave
emission.

This independent (observational) knowledge is important for empirically
constraining QPP models involving the reconnection, for example, for
coalescence of two magnetic loops (Tajima et al., 1987; Kolotkov et al., 2016)
or for magnetic tuning fork model (Takasao & Shibata, 2016), and also flapping
oscillations of the macroscopic current sheet above the loop arcade (Artemyev
& Zimovets, 2012).

In addition, the results obtained are relevant to the interpretation of multi-
modal QPPs in the solar (e.g., Melnikov et al., 2005) and stellar (e.g., Lang
et al., 1983) flares, as the observed second, and sometimes third time
harmonics can correspond to the same monochromatic MHD mode. However, the
difference between the second harmonic obtained in our study from the second
harmonic of the eigen sausage mode is that the non-linearity of the microwave
response leads to its period equals exactly to the half-period of the
principal harmonic, while for the eigen sausage mode it is not so, because of
the wave dispersion (e.g., Inglis & Nakariakov, 2009).

We should stress that our model is deliberately simple, aiming to demonstrate
the discussed effect. A detailed analysis needs to account for the 3D geometry
(see, e.g., Lopin & Nagorny, 2015a; Chen et al., 2015; Lopin & Nagorny, 2015b;
Kolotkov et al., 2021), the beam size of a radio instrument, variation of the
GS spectrum, wave dispersion, etc. An important research avenue is forward
modelling of the manifestation of sausage oscillations in the data delivered
by the future Square Kilometre Array (SKA) instrument, which is expected to
well resolve these oscillations in both time, space and GS spectrum
(Nakariakov et al., 2015).

## Acknowledgements

This study is supported by the Russian Science Foundation project No.
21-12-00195.

## Data Availability

The results obtained in the paper are theoretical; no real observations have
not been used. Distributions of parameters in the magneto-plasma slab were
obtained using equations in Section 2, and their digital version is available
on request to the corresponding author.

## References

  * Antolin & Van Doorsselaere (2013) Antolin P., Van Doorsselaere T., 2013, A&A, 555, A74
  * Artemyev & Zimovets (2012) Artemyev A., Zimovets I., 2012, Sol. Phys., 277, 283
  * Chen et al. (2015) Chen S.-X., Li B., Xiong M., Yu H., Guo M.-Z., 2015, ApJ, 812, 22
  * Cooper et al. (2003) Cooper F. C., Nakariakov V. M., Williams D. R., 2003, A&A, 409, 325
  * Dulk & Marsh (1982) Dulk G. A., Marsh K. A., 1982, ApJ, 259, 350
  * Fleishman & Kuznetsov (2010) Fleishman G. D., Kuznetsov A. A., 2010, ApJ, 721, 1127
  * Gruszecki et al. (2012) Gruszecki M., Nakariakov V. M., Van Doorsselaere T., 2012, A&A, 543, A12
  * Hornsey et al. (2014) Hornsey C., Nakariakov V. M., Fludra A., 2014, A&A, 567, A24
  * Inglis & Nakariakov (2009) Inglis A. R., Nakariakov V. M., 2009, A&A, 493, 259
  * Inglis et al. (2008) Inglis A. R., Nakariakov V. M., Melnikov V. F., 2008, A&A, 487, 1147
  * Kolotkov et al. (2016) Kolotkov D. Y., Nakariakov V. M., Rowlands G., 2016, Phys. Rev. E, 93, 053205
  * Kolotkov et al. (2021) Kolotkov D. Y., Nakariakov V. M., Moss G., Shellard P., 2021, MNRAS, 505, 3505
  * Kupriyanova et al. (2016) Kupriyanova E. G., Kashapova L. K., Reid H. A. S., Myagkova I. N., 2016, Sol. Phys., 291, 3427
  * Kupriyanova et al. (2020) Kupriyanova E., Kolotkov D., Nakariakov V., Kaufman A., 2020, Solar-Terrestrial Physics, 6, 3
  * Kuznetsov & Fleishman (2021) Kuznetsov A. A., Fleishman G. D., 2021, ApJ, 922, 103
  * Kuznetsov et al. (2015) Kuznetsov A. A., Van Doorsselaere T., Reznikova V. E., 2015, Sol. Phys., 290, 1173
  * Lang et al. (1983) Lang K. R., Bookbinder J., Golub L., Davis M. M., 1983, ApJ, 272, L15
  * Li et al. (2020) Li B., Antolin P., Guo M. Z., Kuznetsov A. A., Pascoe D. J., Van Doorsselaere T., Vasheghani Farahani S., 2020, Space Sci. Rev., 216, 136
  * Lopin & Nagorny (2015a) Lopin I., Nagorny I., 2015a, ApJ, 801, 23
  * Lopin & Nagorny (2015b) Lopin I., Nagorny I., 2015b, ApJ, 810, 87
  * Melnikov et al. (2005) Melnikov V. F., Reznikova V. E., Shibasaki K., Nakariakov V. M., 2005, A&A, 439, 727
  * Mossessian & Fleishman (2012) Mossessian G., Fleishman G. D., 2012, ApJ, 748, 140
  * Nakariakov & Melnikov (2009) Nakariakov V. M., Melnikov V. F., 2009, Space Sci. Rev., 149, 119
  * Nakariakov & Roberts (1995) Nakariakov V. M., Roberts B., 1995, Sol. Phys., 159, 399
  * Nakariakov & Verwichte (2005) Nakariakov V. M., Verwichte E., 2005, Living Reviews in Solar Physics, 2, 3
  * Nakariakov et al. (1997) Nakariakov V. M., Roberts B., Petrukhin N. S., 1997, Journal of Plasma Physics, 58, 315
  * Nakariakov et al. (2003) Nakariakov V. M., Melnikov V. F., Reznikova V. E., 2003, A&A, 412, L7
  * Nakariakov et al. (2012) Nakariakov V. M., Hornsey C., Melnikov V. F., 2012, ApJ, 761, 134
  * Nakariakov et al. (2015) Nakariakov V., et al., 2015, in Advancing Astrophysics with the Square Kilometre Array (AASKA14). p. 169 (arXiv:1507.00516)
  * Nakariakov et al. (2018) Nakariakov V. M., Anfinogentov S., Storozhenko A. A., Kurochkin E. A., Bogod V. M., Sharykin I. N., Kaltman T. I., 2018, ApJ, 859, 154
  * Pascoe et al. (2007) Pascoe D. J., Nakariakov V. M., Arber T. D., 2007, Sol. Phys., 246, 165
  * Pascoe et al. (2009) Pascoe D. J., Nakariakov V. M., Arber T. D., Murawski K., 2009, A&A, 494, 1119
  * Reznikova et al. (2015) Reznikova V. E., Van Doorsselaere T., Kuznetsov A. A., 2015, A&A, 575, A47
  * Tajima et al. (1987) Tajima T., Sakai J., Nakajima H., Kosugi T., Brunel F., Kundu M. R., 1987, ApJ, 321, 1031
  * Takasao & Shibata (2016) Takasao S., Shibata K., 2016, ApJ, 823, 150
  * Yu et al. (2013) Yu S., Nakariakov V. M., Selzer L. A., Tan B., Yan Y., 2013, ApJ, 777, 159
  * Zimovets et al. (2013) Zimovets I. V., Kuznetsov S. A., Struminsky A. B., 2013, Astronomy Letters, 39, 267
  * Zimovets et al. (2021) Zimovets I. V., et al., 2021, Space Sci. Rev., 217, 66


# PRoST: Distributed Execution of SPARQL Queries Using Mixed Partitioning
Strategies

Matteo Cossu  Institute for Computer Science
University of Freiburg79110<EMAIL_ADDRESS>,  Michael Färber  Institute for
Computer Science
University of Freiburg79110<EMAIL_ADDRESS>and  Georg
Lausen  Institute for Computer Science
University of Freiburg79110<EMAIL_ADDRESS>
###### Abstract.

The rapidly growing size of RDF graphs in recent years necessitates
distributed storage and parallel processing strategies. To obtain efficient
query processing using computer clusters a wide variety of different
approaches have been proposed. Related to the approach presented in the
current paper are systems built on top of Hadoop HDFS, for example using
Apache Accumulo or using Apache Spark. We present a new RDF store called PRoST
(Partitioned RDF on Spark Tables) based on Apache Spark. PRoST introduces an
innovative strategy that combines the Vertical Partitioning approach with the
Property Table, two preexisting models for storing RDF datasets. We
demonstrate that our proposal outperforms state-of-the-art systems w.r.t. the
runtime for a wide range of query types and without any extensive precomputing
phase.

††copyright: rightsretained††doi: ††isbn: 978-3-89318-078-3††conference: 21st
International Conference on Extending Database Technology (EDBT); March 26-29,
2018; Vienna, Austria††journalyear: 2018

##  1\. Introduction

Organizations and institutions increasingly see a need to represent their data
in a semantically-structured way (Färber et al., rint) and, thus, rely on the
RDF data model. As the data represented in RDF constantly grows in size,
storing and querying these very large RDF graphs becomes a major challenge. In
order to increase the query efficiency, many different approaches have been
proposed e.g. (Abadi et al., 2007; Lee and Liu, 2013; Punnoose et al., 2012;
Schätzle et al., 2016). Most of the current solutions use DBMS-based systems
and map SPARQL queries to SQL queries for the retrieval. Besides the systems
which are implemented for single-node machines, such as RDF-3X (Neumann and
Weikum, 2008), Sesame (Broekstra et al., 2002) and Jena (Wilkinson, 2006),
systems designed for distributed environments are increasingly being used.
They allow to scale up the storage and to provide parallel query execution
capabilities, which is essential to handle very large data. Systems of
distributed environments use in particular Hadoop technologies: S2RDF
(Schätzle et al., 2016) and SPARQLGX (Graux et al., 2016) are implemented on
top of Apache Spark, Rya (Punnoose et al., 2012) on top of Apache Accumulo,
and Sempala (Schätzle et al., 2014) on top of Impala.

Typically, these systems are optimized for specific query patterns and their
data loading times are often sacrificed for better querying performance. Thus,
there is need for a distributed RDF store with better performance on a wide
range of query types, without renouncing a rapid loading phase.

In this work, we describe the approach of PRoST111The source code of PRoST is
online available at https://github.com/tf-dbis-uni-freiburg/PRoST.
(Partitioned RDF on Spark Tables), that aims to improve the efficiency of
queries on RDF data, thereby covering a wide range of query types. Instead of
building a standalone system, such as (Potter et al., 2016), AdPart (Abdelaziz
et al., 2017) or TriAD (Gurajada et al., 2014), PRoST is based on reliable
Hadoop technologies for storing and processing the data and therefore its
performance depends more on the general approach than on the details of the
implementation (e.g. memory management, network protocols). In this paper, we
focus on querying RDF efficiently by means of Apache
Spark222https://spark.apache.org., as it emerged as an important component of
the Hadoop ecosystem: It offers general purpose data processing functions that
exploit efficiently the main memory — differently from MapReduce that relies
mainly on disk operations. This characteristic makes Spark very fast in
practice and able to compute complex queries on large RDF graphs. In
particular, it is possible to obtain high querying performances using Spark
SQL, a module of Spark that adds well-established DBMS techniques to the
distributed computation framework.

Our main contribution is proposing a new storage model for loading RDF graphs
into Hadoop, including an appropriate strategy for translating SPARQL queries
into Spark execution plans. We show that combining the so far separately
applied strategies for representing RDF data in Hadoop — Vertical Partitioning
(Abadi et al., 2007) and the strategy of using Property Tables (Wilkinson,
2006) — can be implemented with relatively little effort and leads to superior
querying performances in many cases, as our evaluation results show.

The paper is structured as follows: In Section 2, we present an overview of
similar systems for querying RDF. In Section 3, we outline our approach for
modeling the data, our strategy for translating SPARQL queries, and the
techniques for improving the performances of PRoST. We summarize our
evaluation results in Section 4 and, in Section 5, we conclude and present
ideas for future implementations of PRoST.

##  2\. Related work

Due to space limitations, in the following we restrict ourselves to Hadoop-
based RDF querying approaches, as they are the most similar to our system
PRoST. S2RDF (Schätzle et al., 2016) is a "SPARQL processor" based on Spark
SQL. It introduces a data partitioning approach that extends Vertical
Partitioning (Abadi et al., 2007) with additional tables, containing
precomputed semi-joins. As a consequence, many intermediate results of queries
are already computed and can be used to shorten the retrieval time. However,
S2RDF trades off the performances with disk space and loading time. For
datasets with a large number of properties (e.g., DBpedia (Auer et al.,
2007)), the time required may make the loading unfeasible.

SPARQLGX (Graux et al., 2016) is another system for distributed SPARQL queries
that uses Vertical Partitioning. SPARQLGX compiles the queries directly into
Spark operations. The system relies entirely on its own statistics to optimize
the computation, in particular for the order of the joins. Differently from
S2RDF and PRoST, SPARQLGX does not use Spark SQL.

Rya (Punnoose et al., 2012) is a popular RDF management system, developed by
an active open-source community and currently an incubating project at Apache.
It is built on top of Apache Accumulo333https://accumulo.apache.org/., a
distributed key-value datastore for Hadoop. Since Accumulo keeps all its
information sorted and indexed by key, Rya stores whole RDF triples as keys.
For this reason, single triples or short ranges can be accessed very fast.
However, considering that the triple table model of Rya contains three
columns, it requires to replicate the data multiple times in order to exploit
the indexes over all the possible elements. Rya uses some query optimization,
e.g. joins reordering, but it lacks of the powerful in-memory data processing
that make, in practice, other systems faster.

##  3\. Approach

The main idea behind PRoST is to store the data twice, partitioned in two
different ways. Each one of the two data structures can be more beneficial (in
terms of performance) than the other regarding certain query types, mainly
because of the different location of the data in the cluster. Our system tries
to exploit this characteristic. We do not limit ourselves to choose the single
best model out of these two for each query, but we split a query in several
parts, such that each sub-query can be executed with the most appropriate
approach. At the end, all the parts can be joined together to produce the
final result. In this way, our approach achieves better performance than when
running only on one of the two representations, as it leverages the synergy
effects of both models.

###  3.1. PRoST Data Model

The two data structures we use to store an RDF graph are Vertical Partitioning
and Property Table. Vertical Partitioning (VP) is an approach proposed by
Abadi et al. (Abadi et al., 2007). It has a strong popularity among RDF
storage systems and it is the main choice of the state-of-the-art approaches
S2RDF (Schätzle et al., 2016) and SPARQLGX (Graux et al., 2016). The main
concept is to create a table for each distinct predicate of the input graph,
containing all the tuples (subject, object) that are connected by that
predicate. The vertical partitioning approach is at the same time powerful and
simple, but it has the drawback that it needs a number of joins proportional
to the number of triple patterns in the query (precisely, the number of triple
patterns minus one). Even with this disadvantage, because the VP tables are
narrow and often small, VP is a valid and practical choice.

The Property Table (PT) is a data scheme that was introduced in the
implementation of Jena2 (Wilkinson, 2006), an influential toolkit for RDF
graph processing, and Sempala (Schätzle et al., 2014), another SPARQL
processor built on top of the Impala framework. A PT consists of a unique
table where each row contains a distinct subject and all the object values for
that subject, stored in columns identified by the property to which they
belong.

The most important advantage of PT, when compared with VP, is that several
joins can be avoided when some of the triple patterns in a query share the
same subject. In this case, part of the query can be executed by a simpler
select operator. Therefore, this approach is beneficial for queries where all
the triple patterns have the same subject variable, usually called star
queries. However, as observed by Abadi et al. (Abadi et al., 2007), the PT has
some issues, in particular the number of NULLs and the multi-valued
attributes. Since not every possible pair subject-predicate has a
corresponding object value, the PT potentially contains a very large number of
NULLs. We solve this problem by storing the table in Parquet
444http://parquet.apache.org, a format that uses run-length encoding. The
other problem of the PT is the presence of multi-valued properties, i.e., when
more than one different object value exists for at least one subject. In this
case, the values are stored using lists that need to be flattened when
executing operations on these attributes. This particular situation produces
an overhead in respect to the VP model, significant if we use the PT for a
single attribute, but negligible if compared to the benefit of avoiding more
than one join to compute the same query.

We have several advantages (e.g. efficient compression) in using the columnar
format Parquet to store the PT, but we then introduce also the undesirable
possibility that the data belonging to a single subject could be split into
different nodes. In order to keep benefiting from the structure of the
property table, we partition horizontally on the subject column. In this way,
we ensure that every row is stored entirely in the same node.

###  3.2. Query Strategies

Having two data structures at the same time requires an abstraction to
represent the queries. For this reason, we define an intermediate format that
we call Join Tree. Each node of this tree represents a sub-query extracted
from a VP table or from the PT. Using Spark, PRoST can execute a query by
calculating the intermediate results from the nodes and then joining them
together in a bottom-up fashion.

We now explain how we translate SPARQL queries into the Join Tree format. For
the sake of clarity, we consider queries with a unique basic graph pattern
without filter, which are a conjunction of triple patterns. The triple
patterns with the same subject are grouped together and translated into a
single node with a special label, denoting that we should use the property
table in this case. All the other groups with a single triple pattern are
translated to nodes that will use the vertical partitioning tables. An example
is shown in Figure 1. There could be a join edge between every pair of nodes
that shares a variable, it means that more than one Join Tree translation is
possible for a single query. Since the tree structure influences the final
performances, choosing carefully the Join Tree is important for the quality of
the system. For this reason, we choose a tree guided by simple statistics, as
explained in the following section.

Figure 1. Example of a conjunction of triple patterns translated to a Join
Tree. The root node uses the Property Table, the others Vertical Partitioning.

###  3.3. Statistics-based Optimization

In the context of our work, joins are the most expensive operators because
they need large portion of the data to be shuffled across the network. The
order of the joins is important to limit this issue, and therefore to speed up
the calculation. An effective way is sorting the joins using statistics of the
input graph. The Join Tree of a query decides indirectly the order in which
the operations are computed. For example, the leaves are computed first, and
the root node is the final one. The statistics we use, simple but effective in
practice, are (1) the total number of triples and (2) the number of distinct
subjects for each predicate. They are calculated during the loading phase
without any significant overhead.

As a general rule of thumb, we want to compute the joins involving few data or
with a high selectivity first (high priority), and finally the ones with many
tuples (low priority). The priority value of a node of the Join Tree is
decided by the following criteria:

  * •

Triple patterns containing literals are scored with the highest priority. The
presence of a literal is a strong constraint that limit the number of
resulting tuples. Therefore, it is a good approach to push down these kind of
nodes.

  * •

A triple pattern of which the underlying data contains many tuples will be
scored proportionally. For instance, a triple with the highest number of
tuples will have also the lowest priority and it will be the root of the tree.
This number is also adjusted according the number of distinct subjects for
that predicate.

  * •

The priority of a node containing data belonging to the Property Table is
scored taking in account all its triple patterns. However, the presence of a
triple pattern with a literal is weighted heavily.

In addition to our effort, Spark SQL’s Catalyst Optimizer (Armbrust et al.,
2015) uses its internal heuristics to improve query performances further. The
trees are not substantially changed, but Spark intervenes in producing
optimized physical plans, since it knows the concrete location of the data on
the cluster. In particular, the optimizer can choose the type of joins to
perform, for example if one of the relations involved is small, a broadcast
join will be performed.

##  4\. Evaluation Results

###  4.1. Benchmark Environment

We perform our tests on a small cluster of 10 machines connected via Gigabit
Ethernet connection. Each machine is equipped with 32GB of memory, 4TB of disk
space and with a 6 Core Intel Xeon E5-2420 processor. The cluster runs
Cloudera CDH 5.11.0 on Ubuntu 14.04 and Spark 2.1.0. Since one machine is the
master, Spark uses only 9 workers. The executor’s memory is 21GB. We use the
WatDiv (Aluç et al., 2014) dataset provided by the Waterloo SPARQL Diversity
Test Suite 7. It is developed in order to measure how an RDF data management
system performs across a wide variety of SPARQL queries with varying
characteristics and has already been applied in the evaluations of the
comparable systems. Given WatDiv we generate a dataset containing around 100
Million RDF triples (5.16 GB) and we evaluate PRoST with the given WatDiv
basic query set. It contains queries of varying shape and selectivity in order
to model different scenarios. The queries are grouped into the following
subsets:

  * •

C: Complex shaped queries.

  * •

F (F1, F2, F3, F4, F5): Snowflake shaped queries.

  * •

L (L1, L2, L3, L4, L5): Linear shaped queries.

  * •

S (S1, S2, S3, S4, S5, S6, S7): Star shaped queries.

###  4.2. Loading

As we can see from Table 1, the storage size used by PRoST decreases when
being stored in HDFS. This is a direct consequence of the compression used by
the Parquet format.

System | Size | Time
---|---|---
PRoST | 2.1 GB | 25m 32s
SPARQLGX | 0.9 GB | 20m 01s
S2RDF | 6.2 GB | 3h 11m 44s
Rya | 3.1 GB | 41m 32s
Table 1.  Size and loading times using WatDiv100M

The database size of S2RDF is the largest in this analysis, since its model
not only needs the creation of VP tables, but also many additional ones
containing the results of its precomputations. SPARQLGX occupies minimal space
because it uses only Vertical Partitioning, where instead PRoST requires
double SPARQLGX’s size for having in addition the Property Table. The loading
time of PRoST is similar to SPARQLGX and around an order of magnitude shorter
than S2RDF. Similarly to what we said before for the database size, S2RDF’s
model requires longer times because of its extensive precomputations.

###  4.3. Vertical Partitioning vs Mixed Strategy

One of our central points is the introduction of the Property Table alongside
Vertical Partitioning. As a consequence, we evaluated the impact of this
addition on query performances.
In Figure 2, we show the times to compare the query set from WatDiv100M with
Vertical Partitioning only and with the additional support of the Property
Table.

Figure 2. Querying time for WatDiv100M with only Vertical Partitioning and
with the mixed strategy. Figure 3. Querying time for WatDiv100M with PRoST,
S2RDF, Rya and SPARQLGX. The scale is logarithmic.

The chart shows clearly that the introduction of the Property Table has a
strong positive impact on performances. For almost every type of query this
version outperforms abundantly the simple Vertical Partitioning approach. The
new strategy is effective for Star queries (S), where most of the triples
share the same subject, as well as Complex (C) and Snowflake (F) queries, that
have more than one different subject. For some of the Linear queries (L), the
results are very similar between the two versions. This result comes from the
fact that Linear queries contain mostly triples with distinct subject
variables, translated using mostly Vertical Partitioning.

###  4.4. PRoST vs Other Systems

To better evaluate our system we compared it in the same environment with
other similar solutions: S2RDF, Rya and SPARQLGX.
In Figure 3 we show the querying times on the WatDiv100M dataset with the
other systems and our implementation. Note that we used a logarithmic scale,
since the large differences would prevent us from visualizing these results
well. Compared to S2RDF, PRoST is faster for the queries F2, S1, S3 and S5 but
slower for the remaining queries. In particular, for the complex queries (C)
and some of the snowflake ones (F3, F4) the S2RDF’s approach makes it faster
by a considerable margin. These differences can be explained with the
extensive precomputations of S2RDF, that heavily decrease the processing time
for joins between VP tables. However we have to keep in mind that S2RDF can
reach these results because of a possibly large loading time, that, for some
datasets, may make it less attractive in practice. In contrast, PRoST relies
on a faster loading phase and its performances does not depend on the
particular input graph, i.e. number of predicates.

Queries | ProST | S2RDF | Rya | SPARQLGX
---|---|---|---|---
Complex | 9,364ms | 3,392ms | 2,195,322ms | 61,363ms
Snowflake | 5,923ms | 1,564ms | 369,016ms | 24,046ms
Linear | 2,419ms | 527ms | 49,044ms | 18,254ms
Star | 1,195ms | 884ms | 6,9606ms | 2,1046ms
Table 2. Average querying time grouped by type of query.

For some queries, Rya is very fast and outperforms PRoST but in other cases it
suffers from significant slowdowns, that are several orders of magnitude
slower. Therefore, its average query time is the worst of the systems
considered, as shown in Table 2. The queries in which Rya performs very well
have in common that their computation involves only few intermediate results.
As a matter of fact, Rya uses powerful indexing methods but it lacks of more
sophisticated approaches to calculate joins. For SPARQLGX, almost all the
queries are computed in around twenty seconds, except the Complex ones (C)
that have a higher number of resulting tuples. PRoST outperforms SPARQLGX in
every case, mostly by around an order of magnitude. Since SPARQLGX uses simple
vertical partitioning, this difference confirms further the improvements
produced by the introduction of the new data organization approach and also
increases the validity of the implementation choices of PRoST. The evaluation
results presented in (Graux et al., 2016) are obtained using a virtual cluster
of 10 nodes on 2 physical machines. In particular, when SPARQLGX is compared
to S2RDF, their results are not consistent with the results from our
experiment, obtained on a cluster of 10 physical machines.

##  5\. Conclusion and Future Work

In this paper, we presented PRoST, a distributed system for RDF storage and
SPARQL querying built on top of Apache Spark. With PRoST we introduce an
innovative data structure for partitioning RDF data that, to the best of our
knowledge, for the first time combines efficiently two pre-existing
approaches, namely Vertical Partitioning and the Property Table. Our
evaluations, in which we compare PRoST on a Hadoop cluster with some state-of-
the-art systems Rya, SPARQLGX, and S2RDF, show that PRoST outperforms Rya and
SPARQLGX in terms of querying time to a large extent and that it achieves
similar results to S2RDF. Notably, PRoST shows consistently good results for
every type of query and is not limited to datasets with certain
characteristics to keep the pre-computation feasible. Therefore, our approach
is in particular appropriate for real-world applications, for which the query
type and the dataset are unknown a priori.

For future work, we considered further improvements in terms of both storage
and querying performance. For example, a promising step might be to add
another Property Table where, instead of the subjects, the rows would be
created around objects. This could be beneficial for triple patterns that
share the same object. Another possible improvement would be to collect more
precise statistics of the input dataset in order to produce better trees and,
hence, a less expensive retrieval.

## References

  * (1)
  * Abadi et al. (2007) Daniel J Abadi et al. 2007\.  Scalable semantic web data management using vertical partitioning. In Proceedings of VLDB. VLDB Endowment, 411–422.
  * Abdelaziz et al. (2017) Ibrahim Abdelaziz et al. 2017\.  Combining Vertex-centric Graph Processing with SPARQL for Large-scale RDF Data Analytics.  IEEE Transactions on Parallel and Distributed Systems (2017).
  * Aluç et al. (2014) Güneş Aluç et al. 2014\.  Diversified stress testing of RDF data management systems. In Proceedings of ISWC. Springer, 197–212.
  * Armbrust et al. (2015) Michael Armbrust et al. 2015\.  Spark SQL: Relational Data Processing in Spark. In Proceedings of ACM SIGMOD. ACM, 1383–1394.
  * Auer et al. (2007) Sören Auer et al. 2007\.  Dbpedia: A nucleus for a web of open data.  The semantic web (2007), 722–735.
  * Broekstra et al. (2002) Jeen Broekstra et al. 2002\.  Sesame: A generic architecture for storing and querying rdf and rdf schema. In Proceedings of ISWC. Springer, 54–68.
  * Färber et al. (rint) Michael Färber et al. 2017 (preprint).  Linked Data Quality of DBpedia, Freebase, OpenCyc, Wikidata, and YAGO.  Semantic Web Journal (2017 (preprint)).
  * Graux et al. (2016) Damien Graux et al. 2016\.  Sparqlgx: Efficient distributed evaluation of sparql with apache spark. In Proceedings of ISWC. Springer, 80–87.
  * Gurajada et al. (2014) Sairam Gurajada et al. 2014\.  TriAD: a distributed shared-nothing RDF engine based on asynchronous message passing. In Proceedings of ACM SIGMOD. ACM, 289–300.
  * Lee and Liu (2013) Kisung Lee and Ling Liu. 2013.  Scaling queries over big RDF graphs with semantic hash partitioning.  Proceedings of VLDB Endowment 6, 14 (2013), 1894–1905.
  * Neumann and Weikum (2008) Thomas Neumann and Gerhard Weikum. 2008.  RDF-3X: a RISC-style engine for RDF.  Proceedings of VLDB Endowment 1, 1 (2008), 647–659.
  * Potter et al. (2016) Anthony Potter et al. 2016\.  Distributed RDF Query Answering with Dynamic Data Exchange. In Proceedings of ISWC. 480–497.
  * Punnoose et al. (2012) Roshan Punnoose et al. 2012\.  Rya: a scalable RDF triple store for the clouds. In Proceedings of 1st IWCI. ACM, 4\.
  * Schätzle et al. (2014) Alexander Schätzle et al. 2014\.  Sempala: interactive SPARQL query processing on hadoop. In Proceedings of ISWC. Springer, 164–179.
  * Schätzle et al. (2016) Alexander Schätzle et al. 2016\.  S2RDF: RDF querying with SPARQL on spark.  Proceedings of VLDB Endowment 9, 10 (2016), 804–815.
  * Wilkinson (2006) Kevin Wilkinson. 2006\.  Jena Property Table Implementation. In Proceedings of SSWKBS.