# EEG2Mel: Reconstructing Sound from Brain Responses to Music

###### Abstract

Information retrieval from brain responses to auditory and visual stimuli has
shown success through classification of song names and image classes presented
to participants while recording EEG signals. Information retrieval in the form
of reconstructing auditory stimuli has also shown some success, but here we
improve on previous methods by reconstructing music stimuli well enough to be
perceived and identified independently. Furthermore, deep learning models were
trained on time-aligned music stimuli spectrum for each corresponding one-
second window of EEG recording, which greatly reduces feature extraction steps
needed when compared to prior studies. The NMED-Tempo and NMED-Hindi datasets
of participants passively listening to full length songs were used to train
and validate Convolutional Neural Network (CNN) regressors. The efficacy of
raw voltage versus power spectrum inputs and linear versus mel spectrogram
outputs were tested, and all inputs and outputs were converted into 2D images.
The quality of reconstructed spectrograms was assessed by training classifiers
which showed 81% accuracy for mel-spectrograms and 72% for linear spectrograms
(10% chance accuracy). Lastly, reconstructions of auditory music stimuli were
discriminated by listeners at an 85% success rate (50% chance) in a two-
alternative match-to-sample task.

##  1 Introduction

Reconstructing stimuli from brain responses has been explored across several
related subdisciplines spanning from signal processing methodology, Cognitive
Neuroscience, and deep learning modeling. A notable example comes from the
culmination of experimental studies showing how to retrieve a Frequency
Following Response (FFR) from averaged brainstem recordings to short ($<$ 1
second) presentations of speech and music [1, 2]. This type of recording only
needs three electrodes, one centered on the scalp, a reference on the earlobe,
and one grounding on the forehead. With this set up, participants are
presented a sound such as a violin playing with their eyes closed repeatedly
for about 1000 - 3000 times. The strength in this approach comes from the
quality of FFRs to be able to contain the fundamental frequency, harmonics,
and onset aligned amplitude from the auditory stimulus. This allows for
researchers to playback an FFR and hear a lower resolution version of the
original stimulus. Such a methodology provided an impressive proof-of-concept
for speech and music retrieval, but it also depends on various experimental
conditions. First, careful stimuli selection is necessary to retrieve FFRs as
it needs to consider stimuli that have clear amplitude bursts, strong onsets,
lower fundamental frequencies ($<$ 300 Hz), and short durations. Furthermore,
the retrieved signal is averaged from hundreds repeated presentations used to
isolate a specific response such as in brainstem activity.

Another exemplary approach in the space of signal processing and Cognitive
Neuroscience comes from auditory tracking experiments where researchers have
been successful in retrieving the speech amplitude from the stimuli [3].
Experiments had participants listening to 180 secs of an audiobook split up
into 25 separate trials. What makes this approach strong is the simplicity
behind the modeling procedure where it takes as input both low frequency and
high-gamma power spectrum signals from Electroencephalogram (EEG) channels and
puts them through a linear model that decodes the input and then maps to the
stimulus target. Despite the modeling simplicity, the approach faces a trade
off with the data processing needing to be more complicated. EEG signals need
to first be split into two signals that separate activity from low and high
frequencies while also being iteratively integrated across channels into the
model over several time lags from 0-250 ms. Nevertheless, this research moves
forward the methodology of neural decoding and acoustic information retrieval
from cortical activity towards a more naturalistic scenario where it does not
require mass averaging of presentations and stimuli presentations are longer.

Figure 1: Visualization of music reconstruction in our study. Brain responses
from music listening are processed by deep regressors and retrieved music is
played back to new participants.

Deep learning methods of decoding stimuli from brain responses have been
successfully implemented through various classification models. In the image
stimuli domain, researchers collected passive cortical responses of people
watching images from ImageNet [4]. They trained the EEG responses in a Long
Short-Term Memory (LSTM) based Recurrent Neural Network (RNN). Model
validation demonstrated that brain responses could be used to correctly
identify an image class at 84% performance that is being perceived (2.5%
chance), i.e., if a participant was watching an image of a panda the model
could correctly identify that the image being perceived was of a panda. In the
auditory domain, recent classification studies have developed methods to
efficiently use EEG responses in Convolutional Neural Networks (CNNs) to
correctly classify the name of the song a participant is listening to. This
was shown to work strongly (85%) when participants listened to music they were
familiar with (8.33% chance), and if the input of the models were formatted as
2D representations of the power spectrum density across channels [5]. Using
different datasets, another study replicated the prior study’s results and
further evaluated the efficacy of EEG data being used as images in computer
vision models implemented through a custom version of AlexNet architecture and
in transfer learning with a pretrained ResNet model. The researchers trained
models with the power spectrum density EEG representations and compared it to
the raw EEG voltage representations. Training on the raw input representation
was able to show that not only was end-to-end classification possible, but
that it could also achieve State-of-the-Art performance at 88.6% with 10%
chance on unfamiliar music [6].

Layer Type | Filter Size | Input | Output
---|---|---|---
Conv2D | $4\times 4$ | $63\times 125\times 1$ | $63\times 125\times 8$
BatchNorm2D | - | - | -
Conv2D | $4\times 4$ | $63\times 125\times 8$ | $63\times 125\times 16$
BatchNorm2D | - | - | -
Conv2D | $4\times 4$ | $63\times 125\times 16$ | $63\times 125\times 32$
BatchNorm2D | - | - | -
Conv2D | $4\times 4$ | $63\times 125\times 32$ | $63\times 125\times 64$
BatchNorm2D | - | - | -
Conv2D | $4\times 4$ | $63\times 125\times 64$ | $32\times 63\times 128$
BatchNorm2D | - | - | -
MaxPool2D | $2\times 2$ | $32\times 63\times 128$ | $17\times 32\times 128$
Flatten | - | $17\times 32\times 128$ | 69632
FC1 | - | 69632 | 128
BatchNorm1D | - | - | -
FC2 | - | 128 | 5632
Reshape | - | 5632 | $44\times 128$

Table 1: Architecture used in our deep regressors. This specific model was
trained on the NMED-H dataset with a spectral input and mel-spectra music
target.
Figure 2: On the left are the two types of input representations we test, and
on the right the two types of target representations for a total of 4 model
combinations as labeled by each arrow. All representations come from
Participant 1 at the 100th second.

Studies moving beyond classification approaches and into trying to reconstruct
the stimuli have seen some initial success with the use of generative models
and careful feature extraction. Following the work in [7], the researchers
transferred weights from their classification models to use as encoders in a
Variational Autoencoder (VAE) and Generative Adversarial Network (GAN). They
successfully recreated the image classes that participants were looking at
with an EEG response as input to the network. For example, if a participant
was looking at images of pandas, then the models would recognize that it is a
panda and output an image of a panda. The novelty of generative models comes
from the ability to not just have outputs that are semantically similar, as
their study claims, but to also generate a distribution of images not directly
presented to participants. This was made possible by the transfer learning of
weights trained to classify image classes, which when fine-tuned to new
examples created training procedures of generative models estimating the
density of image classes. For the purposes of direct information retrieval
this also becomes a limitation because the outputs are estimations from image
classes rather than a pixel-by-pixel reconstruction. This means that if a
participant was looking at a picture of a panda who was sitting down and
profiled to the left, then the model would not necessarily pick up on those
features but return a sampled panda image from the estimated image class
fitted by the model. In music retrieval, another study had found a way to
leverage the acoustic features of the music stimuli and use it as inputs for a
multi-view approach in a deep Variational Canonical Correlation Analysis
(VCCA) model [8]. Despite the lack of quantitative model validation analyses,
this study demonstrated the possibility of using deep learning to be able to
reconstruct music spectra from long EEG responses. In their qualitative
analyses, they confirmed that reconstructed spectra were able to retain
acoustic features from the original stimuli such as pitch, timbre, and tempo.
Furthermore, it allowed for an approach that would move beyond stimuli classes
and attempt to consider time-aligned activity in the stimulus presentation to
be retrieved from brain responses. Here, we present results that advance the
methodology of music retrieval and stimuli reconstruction from brain responses
by training models that have time aligned EEG responses to the music spectra
target. This is done without needing to integrate acoustic features from the
stimuli as input to the model or depending on a multi step feature extraction
process from the EEG. We also present a series of quantitative validation
methods to measure the success of music reconstruction, including feedback
from participants listening to the reconstructed music from our model outputs.
Figure 1 outlines the scope of information retrieval, starting with music
stimuli presented to participants being encoded in their EEG responses and
mapped to the time-aligned stimuli spectrogram by a deep regressor which
outputs a reconstructed music spectrogram that is inverted into a waveform and
presented to a new person.

##  2 Methods

###  2.1 Datasets

Here we train and validate models with the Naturalistic Music
Electroencephalogram Dataset – Tempo (NMED-T) and Naturalistic Music
Electroencephalogram Dataset – Hindi (NMED-H) [9, 10]. Both datasets collect
long recordings of passive brain responses to participants listening to music.
The NMED-T contains recordings from twenty participants listening to ten
different songs that were selected by the researchers because of the lack of
familiarity and differences in tempo. On the other hand, the NMED-H has
recordings from twelve participants listening to four pop songs in Hindi.
These publicly available datasets have been central to various signal
processing studies along with some studies looking into music retrieval from
cortical activity [8, 11]. To aid with replication we provide guiding python
notebooks as examples to our methods and analyses for others to follow. The
supplementary materials also provide examples for others to listen to and
visually inspect. Both datasets provide preprocessed versions of the data
which include standard corrections for faulty channels, line-noise filtering,
and muscle artefact corrections [9, 10]. Here we use those preprocessed
versions of the data, as our focus is on passive cortical signals and not on
any correlated motor behaviors that may show up as muscle artefacts.

###  2.2 Model & Training

Our modeling approach contrasts with previous publications showing success
with decoding and reconstructing complex stimuli, where generative models such
as VAEs and GANs reconstructed image classes and where VCCA was able to
reconstruct music stimuli spectra at the expense of needing a multi-view of
extracted features [7, 8]. We choose a straightforward approach with a
sequential CNN based regressor mapping EEG input directly to the time-aligned
music spectra. Table 1 shows a summary of how the model architecture is
constructed; the models contain five convolutional layers with the last
convolutional layer being the first layer that reduces the dimensionality of
the input. A Max Pooling layer with a small pool size is chosen over a Global
Average Pooling layer used previously during EEG classification [6] using the
same dataset, because it was necessary to limit shrinkage since the output
layer was the size of the spectral target.

Regularization in the model was implemented through dropout layers and
maintaining the size of intermediate fully connected layer as small as
possible. Dropout layers between convolutional layers were kept with 10%
dropout and the dropout layer between the fully connected layers was kept at
15%. The amount and strength of dropout layers were tested demonstrating no
evident effect with less layers and weaker layers, stronger values were not
helpful between convolutional layers but did show some improvement between
fully connected layers. Activation and Kernel L2 regularization was tried but
opted out due to showing signs of a vanishing gradient. Number of filters for
convolutional layers were kept at a base 8 while increasing by a power of 2
for every subsequent layer. Increasing the number of filters showed training
loss outpacing validation loss resulting in overfitted model runs, while
decreasing the number of filters made training loss stagnate too early which
resulted in underfitting. For all intermediate layers a Rectified Linear-Unit
activation was used with a linear output activation. Non-monotonic activation
functions Swish and Mish were tried due to their success in improving image
processing in deep networks, but they did not present any evident advantage
over ReLu during training. Lastly, we use Adaptive Moment Estimation (Adam) as
an optimizer with a 0.0015 learning rate and initialized weights with a He
uniform distribution which have shown advantages in similar training
procedures [6, 12].

The first four minutes of all recordings were used and cut up into five second
chunks. To balance the train and test set distributions across time, we assign
every other chunk to either train or test at a 75/25 ratio. Then all chunks
were split into 1 second examples and randomly shuffled for training and
validation. Five second chunks were also useful in securing consecutive 1
second examples to be reconstructed. These reconstructed five second music
spectra were inverted into waveforms and used as examples in a behavioral
experiment to validate the quality of brain to music reconstruction.
Generalization stays within participants and across unseen chunks of time that
are balanced to be sampled from the beginning, middle, and end of the songs.
Other EEG studies also keep generalization within participants [13, 5], and a
recent study has shown that weak correlations across participants in the
NMED-T could be why generalizing to unseen participants is difficult without
many recorded participants in the dataset [14].

##  3 Experimental Results

###  3.1 Representations

Using the NMED-T, four reconstruction models were trained and evaluated
depending on their varying input and target representations. Figure 2 outlines
the combinations between input and target representations along with their
visual qualities. Each arrow in the figure is a model showing the designated
input to target mapping. Prior classification studies have shown a preference
for input representations either being raw voltage or a Power Spectral Density
(PSD), and their ability to boost performance in CNNs [5, 6]. As a regression
task, we also find it important to compare target representations since linear
and mel-spectrogram representations have tradeoffs. The linear spectrum scales
low frequency activity equally to higher frequency, which allows it to be more
easily inverted back into a listenable waveform. On the other hand, the mel-
spectrogram applies a non-linear scaling across frequency ranges that has made
it easier to map in classification tasks [15, 16], but also adds more noise
during the inversion to a listenable waveform.

Inputs (1 secs) | Targets (1 secs) | Accuracy
---|---|---
PSD (63,125) | Mel-Spec (44,128) | 80.80%
PSD (63,125) | Lin-Spec (44,1025) | 72.28%
Raw (125,125) | Mel-Spec (44,128) | 46.07%
Raw (125,125) | Lin-Spec (44,1025) | 37.53%

Table 2: Summary of reconstruction models’ output classification across the
four representation combinations. Each representation is shown with its data
shape for 1 sec in parentheses.

The reconstruction models were all trained until they reached a plateau in
their loss with varying epoch ranges ( 100-300) needed for the models to
converge their training and validation performance. Given that the mean-
squared error (MSE) loss was only indicative of the convergence of the model
during training and not a metric that informs how image reconstruction
compares across models, we decided to classify the output as an objective
quantitative measure. Simply put, the four trained regressors were used to
create new datasets from their outputs given the original NMED-T as input and
keeping the same train-test split. Then those new output datasets were passed
through a classifier that would attempt to classify the name of the song from
each spectral image reconstruction. This type of validation has been used
before in the image stimuli domain when trying to objectively test the quality
of reconstructed image classes from generative models [7]. The classifiers
were CNN based following the architecture from a previous study where the EEG
was classified into song name classes [6]. Table 2 demonstrates a summary of
the results from the classified outputs across modeling approaches. The raw
input representation showed to not give the best results in our deep CNN
regressors but it did remain well above the chance performance rate of 10%.
Meanwhile, the mel-spectrogram worked better across both the raw and PSD input
representations with performance comparable to studies where only EEG
classification was conducted [17, 13, 18].

Figure 3: Distributions of SSI and PSNR scores across target representations,
along with their mean score.
Figure 4: SSI and PSNR scores across the 128 mel frequency bins of all test
examples. On the right is a reconstructions example zooming in to the lower
and higher frequencies.

Classifying outputs to their stimuli classes reveals fidelity of semantically
relevant reconstructions. This was especially useful in the image stimuli
domain when the models being evaluated focused at the stimuli class level [7].
Furthermore, their results of classification stayed well above chance (2.5%
chance) for forty classes with performance 4̃0%. The results from Table 2 also
stand well above chance for 10 classes (10% chance). Given that our training
method maps EEG to the target stimuli directly rather than estimating stimuli
class densities in the model’s latent space, common image reconstruction
metrics were used to see how different each reconstruction was from the target
image. Structural Similarity Index (SSI) and Peak Signal to Noise Ratio (PSNR)
were chosen because of the ability to compare the reconstruction of perceived
image features and how much noise the reconstructions contain from a lossy
compression comparison respectively [19]. Figure 3 summarizes the results from
both metrics across the test set of the image reconstructions from the
regressor models trained on the PSD input representation along with their
linear and mel-spectrogram targets. SSI measures, as a normalized ratio, the
similarity between images where an SSI score of 1 points to the images being
identical. The mean SSI is higher in the mel-spectrogram reconstruction than
with the linear-spectrogram by 14%. PSNR on the other hand, measures as a
ratio the logged maximum image power over its mean squared error which gives a
relative difference between images and not a normalized score. We also see in
Figure 3 that mel-spectrogram has a higher mean PSNR which further validates
how the mel-spectrogram makes the EEG to music spectra reconstruction closer
to the target.

To better understand what musical features the model learned from mapping to
the mel-spectrogram, we took SSI and PSNR scores across frequencies. Figure 4
demonstrates the results of this analysis; for every example we paired up each
1 second frequency bin from target to model output to calculate both an SSI
and PSNR value. Boxplots were created for all 128 frequency bins, and Figure 4
shows that frequencies vary around a 0.4 SSI score. Meanwhile, the PSNR plot
shows notably lower values for the first 5 frequency bins. We show a visual
example of this difference by taking a reconstruction and zooming into the ten
lowest and highest frequency bins. This visual comparison demonstrates how
many more pixel differences the lower frequencies have when compared to the
higher frequencies. This went against our expectations as we believed it would
be easier to map to lower frequencies because of pattern regularity [20].
Further testing is needed as future studies adopt this methodology but we
speculate that this could be attributed to these model architectures not
explicitly learning temporal dependencies such as in an RNN.

###  3.2 Spectra to Music

The above-mentioned metrics were useful in evaluating the success of the
spectral image reconstructions, but here we are also interested in whether
that allows for the retrieval of perceptually interpretable music
reconstruction. This means that our spectral image reconstructions should be
good enough to then be processed by common out of the box signal processing
libraries that invert spectrograms to soundwaves, and then be identified by
listeners as matching the original stimulus. For a practical experimental
design and for testing the generalization of our methodology, we decided to
train models and produce reconstructed spectral outputs from the NMED-H.
Instead of a total of ten songs like in the NMED-T, the four songs from the
NMED-H made it easier to test reconstructions across time slices and recorded
participants; the training procedure in the deep CNN regressors stayed the
same. Models trained on this dataset converged sooner and provided a steeper
gradient descent during training which could be attributed to the familiarity
of songs from participants [21]. Because the best performing models using the
NMED-T used the PSD and mel-spectrogram as input-target representations, we
focused our training on the NMED-H with that specific representation pairing.
Figure 5 provides a visual summary of the quality of reconstructions from the
NMED-H.

Figure 5: Five second examples of model spectra predictions (left) and their
reconstructions from spectra to sound wave (right). Examples come from the
10th second across randomly selected participants.

As a tradeoff for efficient modeling procedures and feature mapping within the
proposed methodology, processing the outputs to listenable waveforms faces
several steps of lossy transformations. The first comes from the actual deep
learning models themselves where the input to target mapping assumes that the
music signal is hidden/scrambled in cortical activity and attempts to recover
that into a spectral representation. Second, those model outputs must be
denormalized and set to a decibel range. Lastly, the spectrograms are
transformed into waveforms by approximating the Short-Time Fourier Transform
(STFT) magnitude from the mel power spectrum and reconstruction of the phase
is done using Griffin Lim Algorithm (GLA). The mel-spectrogram in comparison
to the linear is more lossy in the last two steps because mel is scaled as a
non-linear function that represents human auditory perception across
frequencies. Nevertheless, our modeling procedure finds it easier to map mel
targets, so we commit to lossy tradeoffs and adjust to produce listenable
examples to present to participants during the behavioral evaluation. During
the denormalization step we use the following transformation:

$\operatorname{Mel\,dB}=(\operatorname{mel}^{T}\times\operatorname{Max\,dB})-\operatorname{Max\,dB}+\operatorname{ref\,dB}$
(1)

In Equation 1, $\operatorname{mel}^{T}$ is the mel-spectrogram we aim to
denormalize set as a matrix transpose and $\operatorname{Max\,dB}$ is set as a
constant to 100 while $\operatorname{ref\,dB}$ is set as a constant within a
song but varies across songs. To avoid presenting jarring noise to
participants, $\operatorname{ref\,dB}$ becomes a free parameter set to [46 dB,
46 dB, 52 dB, 43 dB] for songs 1-4 respectively. The parameter value was set
simply by randomly sampling 5 second examples and listening to how they would
sound with values between 40 dB – 60 dB. The evaluation of which value
produced the least jarring example was decided by the researcher; therefore we
encourage any other studies following this methodology to play around with
their own constraints relevant to their own studies. In Figure 5 it is
demonstrated how well the waveforms come out despite the lossyness of the
inversion process. The reconstructed spectrograms in Figure 5 shows how
information across frequencies are recovered and how that translates to
alignment in amplitude from target to reconstructed waveform.

The produced listenable examples were then presented to participants in a
listening task to evaluate whether their quality was good enough to be
identifiable. The experiment was a two-alternative AB-X task in which a
participant had to pick between sounds A or B that matches the unlabeled sound
X. The target sound X was a 5 second example from one of the four song stimuli
in the NMED-H. Sounds A or B contained a non-corresponding foil and a
corresponding 5 second reconstruction from our deep regressor model that was
taken from spectra to a waveform. Participants were tested on a total of 24
trials plus 4 practice trials at the start of the experiment. The
reconstructed model output produced a total of forty-eight examples which then
half of the examples were used for one version of the experiment and the other
half for a second version. Each version of the experiment controlled for
balancing presentations across time by picking examples that belonged to the
beginning, middle, and end of the songs. Furthermore, the experiment also
contained the same number of examples per song while controlling for each
trial to have reconstructions coming from different participants from the
NMED-H. Trials were all randomized and foils balanced to minimize the
repetition of reconstructions from a specific song or participant. During the
AB-X discrimination task, participants were given a max of 30 minutes to
finish and were allowed to listen to the sound samples as many times as they
liked before making their choices, while also not taking longer than 1 minute
per presentation. This experiment collected responses from a total of 16
participants with 8 participants for each version of the task. On average,
participants had an 85% success rate (50% chance) with max performance of
95.83% and minimum performance of 66.67%. The average performance lines up
with output classification, adding evidence of not just robust image
reconstruction but also the perceived interpretability of listeners when
inverting the spectrum.

##  4 Conclusion

Here we show that it is possible to reconstruct music presented to
participants using their EEG responses as input to a quality that allows
others to correctly identify it when they listen to the music recovered from a
brain response. The EEG to music stimuli mapping approach allows for
reconstructions to preserve time-dependencies, during concatenated examples,
necessary in the perception of music rather than taking an average of the
stimuli class. Furthermore, no heavy feature extraction from the EEG or
stimuli signals are needed in this process and only requires the EEG input to
be transformed into a power spectrum for high quality reconstructions. With
this we have shown that a computer vision approach has allowed for the
processing of EEG responses to long naturalistic music stimuli as well as the
ability to recover information across all frequency components. Finally, our
results present a unique case for computer vision tasks, where what are
normally treated as time-series (EEG and Music) are represented as images for
more efficient training procedures. This leaves a path open for future
investigations focusing on the efficacy of specific computer vision
architectures being tasked to retrieve acoustic information from noisy signals
like EEG.

## References

  * [1] E. Skoe and N. Kraus, “Auditory brainstem response to complex sounds: a tutorial,” _Ear and hearing_ , vol. 31, no. 3, p. 302, 2010.
  * [2] E. B. Coffey, T. Nicol, T. White-Schwoch, B. Chandrasekaran, J. Krizman, E. Skoe, R. J. Zatorre, and N. Kraus, “Evolving perspectives on the sources of the frequency-following response,” _Nature communications_ , vol. 10, no. 1, pp. 1–10, 2019.
  * [3] S. R. Synigal, E. S. Teoh, and E. C. Lalor, “Including measures of high gamma power can improve the decoding of natural speech from eeg,” _Frontiers in human neuroscience_ , p. 130, 2020.
  * [4] C. Spampinato, S. Palazzo, I. Kavasidis, D. Giordano, N. Souly, and M. Shah, “Deep learning human mind for automated visual classification,” pp. 6809–6817, 2017.
  * [5] D. Sonawane, K. P. Miyapuram, B. Rs, and D. J. Lomas, “Guessthemusic: song identification from electroencephalography response,” pp. 154–162, 2021.
  * [6] A. G. Ramirez-Aristizabal, M. K. Ebrahimpour, and C. T. Kello, “Image-based eeg classification of brain responses to song recordings,” _arXiv preprint arXiv:2202.03265_ , 2022.
  * [7] I. Kavasidis, S. Palazzo, C. Spampinato, D. Giordano, and M. Shah, “Brain2image: Converting brain signals into images,” pp. 1809–1817, 2017.
  * [8] A. Ofner and S. Stober, “Shared generative representation of auditory concepts and eeg to reconstruct perceived and imagined music.” pp. 392–399, 2018.
  * [9] S. Losorelli, D. T. Nguyen, J. P. Dmochowski, and B. Kaneshiro, “Nmed-t: A tempo-focused dataset of cortical and behavioral responses to naturalistic music.” pp. 339–346, 2017.
  * [10] B. Kaneshiro, D. T. Nguyen, J. P. Dmochowski, A. M. Norcia, and J. Berger, “Naturalistic music eeg dataset—hindi (nmed-h),” 2016.
  * [11] A. Vinay, A. Lerch, and G. Leslie, “Mind the beat: detecting audio onsets from eeg recordings of music listening,” pp. 231–235, 2021.
  * [12] M. Ebrahimpour, T. Shea, A. Danielescu, D. Noelle, and C. Kello, “End-to-end auditory object recognition via inception nucleus,” pp. 146–150, 2020.
  * [13] M.-A. Moinnereau, T. Brienne, S. Brodeur, J. Rouat, K. Whittingstall, and E. Plourde, “Classification of auditory stimuli from eeg signals with a regulated recurrent neural network reservoir,” _arXiv preprint arXiv:1804.10322_ , 2018.
  * [14] P. Pandey, G. Sharma, K. P. Miyapuram, R. Subramanian, and D. Lomas, “Music identification using brain responses to initial snippets,” in _ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2022, pp. 1246–1250.
  * [15] S. Gururani, C. Summers, and A. Lerch, “Instrument activity detection in polyphonic music using deep neural networks.” pp. 569–576, 2018.
  * [16] S. Chillara, A. Kavitha, S. A. Neginhal, S. Haldia, and K. Vidyullatha, “Music genre classification using machine learning algorithms: a comparison,” _Int. Res. J. Eng. Technol.(IRJET)_ , vol. 6, no. 05, pp. 851–858, 2019.
  * [17] S. Stober, D. J. Cameron, and J. A. Grahn, “Using convolutional neural networks to recognize rhythm. stimuli from electroencephalography recordings,” _Advances in neural information processing systems_ , vol. 27, 2014.
  * [18] Y. Yu, S. Beuret, D. Zeng, and K. Oyama, “Deep learning of human perception in audio event classification,” in _2018 IEEE International Symposium on Multimedia (ISM)_. IEEE, 2018, pp. 188–189.
  * [19] A. Hore and D. Ziou, “Image quality metrics: Psnr vs. ssim,” pp. 2366–2369, 2010\.
  * [20] T. Lenc, P. E. Keller, M. Varlet, and S. Nozaradan, “Neural tracking of the musical beat is enhanced by low-frequency sounds,” _Proceedings of the National Academy of Sciences_ , vol. 115, no. 32, pp. 8221–8226, 2018.
  * [21] S. K. Hadjidimitriou and L. J. Hadjileontiadis, “Eeg-based classification of music appraisal responses using time-frequency analysis and familiarity ratings,” _IEEE Transactions on Affective Computing_ , vol. 4, no. 2, pp. 161–172, 2013.


# Differentiating-based Vectorization for Sparse Kernels

Zachary Cetinic  University of Toronto
Toronto, Canada
<EMAIL_ADDRESS>Kazem Cheshmi  University of Toronto
Toronto, Canada
<EMAIL_ADDRESS>Maryam Mehri Dehnavi  University of Toronto
Toronto, Canada
<EMAIL_ADDRESS>
###### Abstract

Sparse computations frequently appear in scientific simulations and the
performance of these simulations rely heavily on the optimization of the
sparse codes. The compact data structures and irregular computation patterns
in sparse matrix computations introduce challenges to vectorizing these codes.
Available approaches primarily vectorize regular regions of computations in
the sparse code. They also reorganize data and computations, at a cost, to
increase the number of regular regions. In this work, we propose a novel
polyhedral model, called the partially strided codelets (PSC), that enables
the vectorization of computation regions with irregular data access patterns.
PSCs also improve data locality in sparse computation. Our DDF inspector-
executor framework efficiently mines the memory accesses in the sparse
computation, using an access function differentiation approach, to find PSC
codelets. It generates vectorized code for the sparse matrix multiplication
kernel (SpMV), a kernel with parallel outer loops, and for kernels with
carried dependence, specifically the sparse triangular solver (SpTRSV). We
demonstrate the performance of the DDF-generated code on a set of 60 large and
small matrices (0.05-130M nonzeros). DDF outperforms the highly specialized
library MKL with an average speedup of 1.93 and 4.5$\times$ for SpMV and
SpTRSV, respectively. For the same matrices, DDF outperforms the state-of-the-
art inspector-executor framework Sympiler [1] for the SpTRSV kernel by up to
11$\times$ and the work by Augustine et. al [2] for the SpMV kernel by up to
12$\times$.

###### Index Terms:

Vectorization, Sparse Matrix Computations, Polyhedral Analysis

##  I Introduction

Sparse matrix computations are important kernels used in a large class of
scientific and machine learning applications. The performance of sparse
kernels is noticeably improved if the code is vectorized to exploit single
instruction multiple data (SIMD) capabilities of the underlying architecture.
Vectorization potentially increases opportunities to optimize for locality,
further increasing the performance of the sparse code. SIMD instructions can
efficiently vectorize groups of operations that access consecutive data, i.e.
have a regular access pattern. A computation is regular when its operations
access memory addresses with a constant distance from each other, i.e. strided
access, and otherwise irregular. However, because a sparse matrix is typically
stored in a compact representation [3], the sparse matrix code accesses are
often irregular, making vectorization challenging. Not vectorizing irregular
segments in the sparse matrix calculation can potentially reduce locality and
the performance of the sparse matrix kernel.

The most common approach to improve vectorization of irregular regions in
sparse matrix codes is to reorganize data or computations so the region
becomes regular and amenable to vectorization. These work typically optimize
sparse kernels with no loop-carried dependencies, for example [4, 5, 6, 7, 8,
9, 10] improve the performance of sparse matrix-vector multiplication (SpMV)
with reorganization. CSR5 [11] and CVR [12] reorganize data to create
contiguous accesses and reorganize computations to increase the number of
independent operations for vectorization. However, their computation
reorganization typically increases the number of floating-point operations in
the kernel, reducing the benefits of vectorization. Sparse storage formats
such as ELL and DIA [13] are used to pad the sparse matrix with additional
nonzeros. Padding increases regular regions in the matrix and hence increases
the number of consecutive vector loads. In sparse matrices with very irregular
patterns, the amount of padding will increase, leading to more floating-point
operations. Reorganization approaches are often library-based and thus are
specialized for a specific set of matrix patterns or the access patterns
appearing in a specific matrix. They also need to be manually ported to new
architectures. For example, work in [14] primarily optimizes matrices from
power-law graph computations, CVR is for sparse graphs, and ELL is most
efficient for block-structured matrices such as factors of a direct solver
[15].

Compilers, such as those based on the polyhedral model, generate vectorizable
and portable code for affine codes such as dense matrix computations [16, 17,
18, 19, 20, 21]. However, compilers are limited in vectorizing sparse codes
because of the existing indirection in the memory access patterns of sparse
computations, caused by data compaction. This indirection introduces
additional challenges to compilers when optimizing sparse kernels with loop
carried dependencies, e.g. the sparse triangular solver. Recently the
polyhedral method was extended to optimize sparse matrix computations via
inspector-executor approaches [22, 23, 24]. The inspectors execute at runtime
to resolve indirection and find data dependencies between operations and this
information is used to transform the original source code into optimized code
called the executor.

Figure 1:  Various generated codes for the given sparsity pattern (shown in
the matrix) used to compute the Sparse Matrix-Vector multiplication (SpMV)
kernel. a) Shows the code that computes SpMV using the compressed sparse row
(CSR) storage format. b) Shows the DDF generated code. c) Shows the code
generated from the Sympiler framework [1]. d) Shows the code generated using
the work by Augustine et al.[2]

Amongst inspector-executor frameworks, Augustine et al. [2], which we refer to
as the PIC framework, and Sympiler [1] inspect the sparsity pattern of input
matrices to generate vectorized code for sparse matrix computations. The work
in [2] proposes an inspector to mine the memory trace of the SpMV kernel to
find regular regions in the computation for vectorization. The remaining
regions are unrolled. Unrolling irregular regions limits the scalability of
this work because of increasing code size for large matrices and for matrices
with high irregularity; their framework can primarily optimize small sparse
matrices (with less than 0.5 million nonzeros). Also, they do not support
sparse kernels with loop-carried dependencies. Sympiler uses an inspector to
resolve data dependencies in sparse kernels such as SPTRV. Its inspector finds
iterations that operate on rows with similar patterns, i.e. regular regions,
and applies code transformations that primarily optimize for tiling and
vectorization. Sympiler adds additional nonzeros to some of the rows with non-
consecutive memory accesses and thus pads the input matrix to form row-blocks.
Row-blocks are regular regions that can be easily vectorized. To limit the
overheads from padding, Sympiler only creates row-blocks when profitable, thus
limiting the size of row blocks it creates. It is thus not able to efficiently
optimize sparse matrices with few rows of similar sparsity patterns.

The standard approach of prior work in vectorizing sparse kernels is to either
find regions with regular computations or create more regular regions in the
computations with for example data reorganization. The focus of this work is
to expose vectorization opportunities in computation regions with irregular
computation patterns. We present a novel polyhedral model which we call the
partially strided codelet that increases opportunities for vectorization of
sparse matrix codes and improves data locality. PSCs represent computation
regions that have non-consecutive, i.e. unstrided, memory accesses but can
still benefit from vectorization because a subset of their accesses is
strided. We also present a novel inspector-executor framework called DDF
(Differentiating Data access Functions) that efficiently inspects the memory
access patterns in sparse matrix codes to mine for PSCs using an access
function differentiation approach and to generate optimized vectorized code
for the matrix computation. DDF’s inspector supports sparse kernels with and
without loop-carried dependencies as it is able to look for data dependencies
in the computation region.

The contributions of the work are:

  * •

A novel polyhedral model, called the partially strided codelet, that enables
the vectorization of computation regions with non-consecutive memory accesses
in sparse codes while improving temporal and spatial locality.

  * •

A memory access differentiation approach, that uses the first order partial
difference (FOPD) of access functions to distinguish strided access from non-
strided accesses in the sparse matrix computation. FOPD’s are used to mine for
PSCs in spare matrix codes.

  * •

The DDF inspector-executor framework that mines the data accesses of sparse
codes with or without dependencies to find PSCs. The inspector executes in
parallel with low overhead and generates a compact code that can execute in
parallel.

  * •

We demonstrate the performance of DDF for SpTRSV and SpMV kernels. The speedup
of the DDF-specialized code for SpTRSV is on average 1.79$\times$ and
4.5$\times$ faster than that of Sympiler and MKL respectively. For SpMV, DDF-
generated code is on average 1.42$\times$ and 1.93$\times$ faster than that of
CSR5 and MKL respectively. For small problems that the work from [2] supports,
DDF is 10.57$\times$ faster than their approach.

##  II Motivation

In this section, we use the matrix in Figure 1 to explain the approach used in
inspector-executor frameworks, i.e. Sympiler and the work in [2] (which we
call PIC), to optimize different computation regions in the SpMV kernel and
compare to the approach proposed in this work. The example uses the compressed
sparse row (CSR) format to store matrix $A$ and then computes $y=A\times x$
where $x$ is a vector. The region with the green border is a regular
computation region in the SpMV computation and the red and orange colors show
regions with non-strided accesses.

All tools can efficiently vectorize the green region since all of its accesses
are consecutive and thus strided. BLAS routines are typically used to optimize
such regions. They enable vectorization, and in part, also improve locality by
reusing the entries in the input vector (e.g. x[6], x[7]) and the output
vector (e.g. y[3], y[4]).

For the non-strided segments, Sympiler uses a padding strategy, when
profitable, to create regular regions that can be easily vectorized. However,
Sympiler will not pad the red and yellow regions, since they are too
irregular, and padding will introduce a large number of additional operations.
Instead, Sympiler will fall back to the original non-affine code and miss any
additional vectorization opportunities. PIC will generate fully unrolled code
for the non-strided regions, as shown in Figure 1d. While this reduces the
number of instructions, it increases code size and as a result, increases the
number of instruction cache misses, especially for large matrices, limiting
PIC’s scalability.

For irregular regions, DDF uses partially strided codelets to improve locality
and to provide a better vectorization efficiency. The SpMV code for the red
region has non-strided accesses in x, while its accesses to Ax and y are
strided. DDF generates the PSC codelet bordered with red in Figure 1b. The
code stores the non-strided addresses in s to be reused across all three
iterations of $i$ instead of loading indices for every operation. This also
enables the reuse of index Ap[i] in Ax so it does not have to be loaded per
each outermost iteration. Similarly, for the yellow region, accesses to both
Ax and x are non-strided, but vectorizing computations of the innermost loop
improves spatial locality due to contiguous accesses in Ax.

##  III Partially Strided Codelet

This section introduces partially strided codelets and discusses their
classification as well as efficiency in optimizing computation regions with
irregular access patterns. We also present a novel cost model and a
differentiation-based PSC detection strategy, both of which are used in the
DDF framework to efficiently find PSCs in sparse matrix computations.

Figure 2: Polyhedral representation of the SpMV kernel shown in Figure 1a

###  III-A Definitions

Polyhedral model and data access functions. A loop nest that contains a set of
statements is represented with a polyhedral model through an integer
polyhedron sets $\mathcal{I}$ and relations $f$. A statement is made of a data
space described with $\mathcal{D}$ which is a disjoint set containing
$\mathcal{D}_{0},\ldots,\mathcal{D}_{n}$. An integer polyhedral set
$\mathcal{I}=[i_{0},...,i_{n}]$ is a collection of inequalities that create
bounds for each dimension inside $i\in\mathcal{I}$. For each
$\mathcal{D}_{d}\in\mathcal{D}$ a data access function $f$ is used to describe
how the data space $\mathcal{D}_{d}$ is accessed by the iteration space of
$\mathcal{I}$. In other words, a data access function maps an iteration space
to a data space, i.e. $f_{\mathcal{I}\rightarrow\mathcal{D}}$ . The SpMV code
in Figure 1a has one statement. That statement has three data spaces y, Ax,
and x as well as three data access functions. Figure 2 shows the polyhedron
sets for $\mathcal{I}=[i_{0},i_{1}]$ and also the access functions
corresponding to each data space.

Figure 3: a) Shows the set of SpMV operations from the red outlined region of
the matrix in Figure 1. b) Shows the iteration space and access functions of a
codelet that maps to the data-space seen in section $a)$ of the figure.

Codelet. A polyhedral model that has a convex integer polyhedron with no flow
dependencies, i.e. read after write dependency between access functions, is a
codelet. A codelet only has one statement and the operation in that statement
is a SIMD-supported operation. A set of SpMV operations is shown in Figure 3a
and its polyhedral representation including iteration space and data access
functions are shown in Figure 3b. All operations in this model are independent
and have only one statement (y[*] += Ax[*] * x[*]), with a multiply-add
operation which is supported by Intel and AMD SIMD units. These properties
satisfy the criteria for being a codelet.

Strided data access function. A function that can be expressed with a linear
combination of induction variables in $\mathcal{I}$ is a strided access
function. The data access function $f$ and $g$ in Figure 3b are both strided
and can be expressed as a linear combination of $[i_{0},i_{1}]$.

###  III-B Partially Strided Codelet Classification

An efficient way to vectorize a codelet is to find operations with strided
accesses across different iterations. This enables the vectorization of more
operations and potentially increases data reuse between different iterations.
However, current approaches are primarily limited to vectorizing codelets that
all of their access functions are strided. An efficient BLAS [25]
implementation is used for this purpose and thus, we call these codelets types
BLAS codelets. In this work we define a novel set of codelets called Partially
Strided Codelets that have at most $n-1$ strided access functions and at least
one non-strided access function. These codelets can benefit from vectorization
because they have one or more strided access functions.

PSCs are classified into different types based on the number of access
functions that are strided. For example in SpMV and SpTRSV that have codelets
with three access functions, $n=3$, two types of PSCs can be defined. The PSC
type I codelet is used when two of access functions are strided, and the PSC
type II codelet is used when only one access function is strided. Figure 1b
shows two types of PSCs, outlined in yellow and red, and also the BLAS
codelet, outlined in green.

###  III-C PSC Efficiency and Cost Model

Figure 4: Shows best performing codelet type on a wide variety (150K total
points) of computational regions found in SpMV/SpTRSV.

Depending on the number of operations and strided accesses in a computation
region, either PSC or BLAS codelets can be profitable. To illustrate this, we
extracted over 150,000 computation regions, with a different number of
operations and strided accesses, from the SpMV and SpTRSV kernels on 60
matrices in the Suitesparse repository [26]. Each region is optimized using
one or several codelets of the same type and the best performing codelet type
for that region is inserted as a point in Figure 4. The red, green, and blue
colors are used for PSC I, PSC II, or BLAS codelets respectively. The figure
shows that some regions perform best with PSCs and some with BLAS codelets,
e.g. PSC I provides the fastest execution time over all codelets for regions
with a large percentage of non-strided accesses.

The performance of codelets differs based on the number of memory accesses
that occur from running the codelets. Per Section III-A, an access function is
described as $f(\mathcal{I})=x+s[i_{0},\ldots,i_{n}]$ where $x$ is an integer
and $s$ is an $n$-dimensional array. If $f$ is a strided access function, an
integer $s[i]$ for each $i\in\mathcal{I}$ needs to be loaded from memory
before the operation in the codelets executes. Each integer, $s[i]$, is the
$i^{th}$ coefficient in the linear combination of $f$. However, for a non-
strided access function, the number of memory accesses in the codelet is equal
to the size of $s[i]$ for $i\in\mathcal{I}$, which is larger than one. Thus,
compared to BLAS codelets, PSCs access more memory addresses. But this does
not imply that PSC codelets are inefficient because the number of codelets
used to cover a computational region also affects the total memory accesses.
For example, the computation region in Figure 3b is a PSC I with one non-
strided access function $h(i_{0},i_{1})=i_{0}+s[i_{1}]$ where $s=\\{0,2,5\\}$
and thus requires loading 4 memory addresses. If BLAS is selected for the same
region, a minimum of 3 codelets will need to be used leading to 9 loads.

Figure 5: Correlation between predicted cost of codelets and average execution
time of each codelet. Averages were obtained using 1,000,000 executions per
codelet.

To efficiently determine the types of codelets for a computation region, we
first define a cost model to compute the cost of a single codelet. The cost
model uses the number of memory accesses of each memory access function, i.e.
$|s_{d}|+1$ and the number of operations in a codelet($|p|$) to estimate the
execution time of the codelet. The cost of a codelet $l$ is $c_{l}$ defined
as:

$c_{l}=|p|+m+\sum_{d=0}^{m}|s_{d}|$ (1)

If a computation region is described with $k$ codelet types then the overall
cost of the codelets for that region is $C=\sum_{l=0}^{k}c_{l}$. Figure 5
shows the correlation between the codelet cost model and the execution time of
codelets. We use a set of codelets with different sizes from SpMV and SpTRSV
on different computation regions of the sparse matrices in Suitesparse [27].
The X-axis shows the cost of each codelet, and the Y-axis shows their
corresponding execution time in seconds. As shown, the cost model predicts the
performance of the codelets with a correlation of $(r^{2}=0.99)$.

###  III-D Differentiation Based PSC Detection

In this section, we explain how the first order partial difference (FOPD) of
the access functions in a computation region can be used to detect a codelet
type.

First Order Partial Differentiation (FOPDi). Given the data access function
$f$ with iteration space of $\mathcal{I}=[i_{0},i_{1}]$, the first order
partial difference of $f$ with respect to $i_{1}\in\mathcal{I}$ is computed as
$\frac{\Delta}{\Delta
i_{1}}f(\mathcal{I})=\Delta_{i_{1}}f=f(i_{0},i_{1}+1)-f(i_{0},i_{1})$. FOPD
shows if accesses to a data space are strided with respect to the induction
variable $i_{1}$. Figure 6 illustrates the process of computing the FOPD for
the access function $h$ given the computation region shown in Figure 3. For
example, the FOPD of $h$ evaluated at $i_{0}=1,i_{1}=1$ is
$\Delta_{i_{1}}h(1,1)=3$.

Figure 6: Illustrates the process of taking the derivative of the access
function $h$ with respect to $i_{0}$ and $i_{1}$.

FOPD of access functions are used to distinguish types of codelets by finding
strided access functions. Given a codelet with three access functions and the
iteration space of $\mathcal{I}=[i_{0},i_{1}]$, an access function $f$ is
strided if its FOPDs with respect to $\mathcal{I}$ are equal in the entire
iteration space. In other words, $f$ is strided if the elements in
$\Delta_{i_{0}}f(i_{0},i_{1})$ are equal to each other and similarly for
elements in $\Delta_{i_{1}}f(i_{0},i_{1})$. With the strided definition above,
all codelet types can be defined per definitions in Section III-B. For
example, the first three accesses in function $g(i_{0},i_{1})$ in Figure 6 are
to consecutive locations (0, 1, 2) in Ax. The FOPD of the first two accesses,
wrt. $i_{1}$, is FOPD${}_{i_{1}}(0,0):g(0,1)-g(0,0)=1$ and for the second and
third accesses is FOPD${}_{i_{1}}(0,1):g(0,2)-g(0,1)=1$. Similarly the FOPD of
the accesses for $i_{0}$ are FOPD${}_{i_{0}}(0,0):g(1,0)-g(0,0)=3$ and
FOPD${}_{i_{0}}(1,0):g(2,0)-g(1,0)=3$. Since FOPD${}_{i_{1}}(0,0)$ and
FOPD${}_{i_{1}}(0,1)$ are equal, and FOPD${}_{i_{0}}(0,0)$ and
FOPD${}_{i_{0}}(1,0)$ are equal, the function in the iteration space
$\mathcal{I}=\\{i=0\land 0\leq j\leq 3\\}$ is strided. Because function $f$ is
also strided and $h$ is not strided, the codelet is categorized as a PSC type
I.

⬇

1#include "DDF.h"

2int main(){

3 Kernel SpMV;

4 SpMV.generate_c("spmv.h", Arch::x86);

5 return 1;

6}

Listing 1: The input code to DDF.

⬇

1#include "DDF.h"

2#include "spmv.h"

3int main(){

4 Matrix A("input.mtx");

5 Vector x("x.mtx"), y("y.mtx");

6 /// \------ Inspector \------ ///

7 vector<codelet *> clist=DDF_Inspector(A.pattern(), SpMV);

8 /// \------ Executor \------ ///

9 SpMV_Vectorized(clist, A, x, y);

10 return 1;

11}

Listing 2: A sample driver code to call the DDF-generated code.

##  IV DDF: Differentiating Data Access Functions to Mine PSC Codelets

DDF is an inspector-executor framework that finds partially strided codelets
from an input code and the matrix sparsity pattern using an inspector and
generates a vectorized code as its executor. The DDF input specification is
shown in Listing 1. It takes the input kernel, SpMV or SpTRSV, and the target
architecture type and generates a driver code shown in Listing 2 and a
vectorized code in spmv.h. The vectorized code is generated when generate_c is
called in line 4 of Listing 1. In lines 4–5 of the Listing 2, the input matrix
$A$ and vectors $x$ and $y$ are loaded. DDF’s inspector then creates a list of
codelets in line 7, and the vectorized code executes in line 9.

###  IV-A The DDF Inspector

DDF’s inspector uses Algorithm 1 to generate efficient code for the input
sparse kernel. The algorithm first creates groups of independent iterations
and then for the computation regions in each group finds a best codelet
combination to generate.

#### Inputs and output

The inputs to the DDF inspector algorithm are a kernel code $\mathcal{K}$,
i.e., SpMV or SpTRSV, and the pattern of the input matrix $\mathcal{P}$. The
algorithm generates a list of codelets $cList$ to be used in the vectorized
code. In DDF, the kernel code $\mathcal{K}$ is represented with an abstract
syntax tree (AST) that represents loops and operations. The pattern
$\mathcal{P}$ is stored in a compressed sparse row (CSR) storage.

Input : $\mathcal{K},\mathcal{P}$

Output : $cList$

/* 1) Grouping independent memory accesses */

1 $f_{0},f_{1},f_{2},\mathcal{I}\leftarrow
computeAccessFunctions(\mathcal{K},\mathcal{P})$;

2 $\mathcal{G}\leftarrow findDependencies(f_{0},f_{1},f_{2},\mathcal{I})$;

3 $\mathcal{S}\leftarrow partitionIterationSpace(\mathcal{G},\mathcal{I})$;

/* 2) Codelet Creation */

4 for  _$p\in\mathcal{S}$_ do

5 FOPDs $\leftarrow ComputeFOPD(p,f_{0},f_{1},f_{2})$;

6 $\mathcal{R}\leftarrow GetConsecutiveIterations(p)$;

7 for  _$r\in\mathcal{R}$_ do

8 $clb\leftarrow BLAS\\_first(r,$ FOPDs);

9 $clp1\leftarrow PSCI\\_first(r,$ FOPDs);

10 $clp2\leftarrow PSCII\\_first(r,$ FOPDs);

11 $min\\_cl\leftarrow min\\_cost(clb,clp1,clp2)$;

12 $cList.append(min\\_cl)$;

13

14

Algorithm 1 DDF Inspector

####  IV-A1 Grouping independent memory accesses

The first step of the inspector algorithm computes memory access functions of
the input code and then generates a graph $\mathcal{G}$ and uses it to create
groups of operations that are independent to support parallelism. Function
$computeAccessFunctions$ in line 1 of Algorithm 1 uses kernel $\mathcal{K}$
and goes over $\mathcal{P}$ to compute access functions and the iteration
space of the input kernel. Then $findDependencies$ in line 1 uses access
functions to compute a graph $\mathcal{G}$. The vertices in this graph
represent iterations of the outermost loop $i\in\mathcal{I}=[i,j]$ and its
edges are the flow dependencies across iterations of $i$. In kernels with loop
carried dependencies, $\mathcal{G}$ is passed to the LBC algorithm [28] to
create independent parallel workloads. For sparse kernels with parallel loops,
$\mathcal{G}$ does not have any edges, thus instead, the function groups
iterations based on the number of operations. Finally, iterations of the
outermost loop are grouped and the groups are stored in $\mathcal{S}$ using
$partitionIterationSpace$ in line 3.

####  IV-A2 Codelet creation

For every group of iterations $p\in\mathcal{S}$, the algorithm generates the
FOPD of the access functions in line 5. Function $GetConsecutiveIterations$ in
line 1 groups operations across $t$ consecutive iterations and creates a
number of computation regions stored in $\mathcal{R}$. To reduce the overhead
of mining, the algorithm limits the search window in line 6 to $t$ consecutive
iterations. Our experiments show that this does not lead to slowdowns because
consecutive iterations in sparse codes typically operate on consecutive rows
of a matrix, and thus their vectorization can potentially improve spatial
locality.

Figure 7: Shows generic codes used to evaluate different codelet types; all
codelets are parameterized by m, n which bound the outer and inner dimensions
respectively. Each type also uses f,g and h which are special structs
containing offset information used to process each codelet.

Lines 8-12 in the algorithm inspect each computation region $r$ and finds the
best codelet combination for that region using three different strategies,
i.e. $BLAS\\_first$ , $PSCI\\_first$, and $PSCII\\_first$. Since the
strategies mine for codelets, they take as input the FOPDs computed from line
5 along with $r$. The BLAS-first strategy prioritizes finding BLAS codelets
and then mines for PSCs. The PSC I-first strategy mines for PSC type I
codelets first, and the PSC II-strategy looks for PSC type II codelets and
then other codelets. To reduce the execution time of each strategy, we store
the list of codelets for each mined region and look up this list when
executing the other strategies. The costs of codelets from each approach are
stored in constants $clb$, $clp1$, and $clp2$ and then in line 11 the
combination with the lowest cost is chosen. Line 12 appends the most efficient
codelet combination found for the current computation region to the output
list of codelets in the Algorithm, i.e. $cList$.

###  IV-B DDF Executor

The DDF driver runs the executor code, as shown in line 9 of listing 2. To
efficiently use the instruction cache and to eliminate the need for
recomputing a new executor code per matrix pattern, DDF generates a parametric
code. The parametric code for the three types of codelets, BLAS, PSC I, and
PSC II are shown in Figure 7. As shown each generic code of a codelet class
takes the data access functions as input and can vectorize any codelet of the
same class. These generic codelets are called inside a switch statement in the
executor code, shown in Listing 3. In the executor, each mined codelet in the
codelet list is mapped to a generic codelet using this switch statement. A
change in the input sparsity pattern results in a different codelet list,
however, the same parametric code is still used. For matrices of different
size, the codelets list size changes, however, the same executor code in
Listing 3 can be used. Thus, the code size is invariant to the size of matrix.

⬇

1#include "Codeletsx86.h"

2

3void vectorized_code(vector<Codlet*> clsit, CSR* A, double *x, double *y){

4 for(int i = 0; i < clist.partitions(); i++){

5#pragma omp parallel

6 for(int j = 0; j < clist[i].size(); j++){

7 switch(clist[i][j].type){

8 case BLAS:

9 BLAS_codelet();

10 break;

11 case PSCI:

12 PSCI_codelet();

13 break;

14 case PSCII:

15 PSCII_codelet();

16 break;

17}}}}

Listing 3: The vectorized code for SpMV generated by DDF in spmv.h.

##  V Results

We evaluate the performance of DDF using two kernels, sparse triangular solver
(SpTRSV) and sparse matrix-vector multiplication (SpMV). DDF is compared to
two other inspector-executor approaches, i.e. Sympiler and PIC [2]. We also
compare DDF to libraries CSR5 [11], MKL [29], and our implementation of SpMV
using the ELLPACK [13] storage format (referred to as ELL in the figures).
Sympiler does not support SpMV, and SpTRV is not supported by the CSR5, PIC,
and ELL implementations. Hence those tools are omitted from the respective
figures.

The set of symmetric positive definite (SPD) matrices is used for evaluation.
These matrices are selected from the SPD symmetric matrices in the SuiteSparse
repository [26] and are diverse in both size and sparsity pattern. For the
SpTRSV evaluations, we only use the lower triangular half of the SPD symmetric
matrices. Matrices with ID 0=30 are the thirty largest matrices in the
repository, Matrix IDs 31-49 are L-factors [30] of the largest matrices, and
10 small matrices (ID = 50-60) with 50-350K nonzeros are also included.
L-factors matrices are included to compare the performance of DDF to Sympiler
because Sympiler is best suited to optimize these matrices. Small matrices are
included to compare with PIC as it does not scale to larger problems. The
test-bed architectures are an Intel(R) Xeon(R) Gold 5115 CPU (2.8GHz, 14080K
L3 Cache) with 20 cores and 64GB of main memory (Intel) and an AMD Ryzen 3900x
CPU (3.8GHz, 64MB L3 Cache) with 12 cores and 32GB of main memory (AMD). The
performance of kernels is reported on both architectures but the analysis is
only shown for the Intel machine. All generated code, implementations of
different approaches, and library drivers are compiled with GCC v.7.2.0
compiler and with the -O3 flag. All benchmarks are executed 5 times and the
median value of all runs is reported. To report the performance in GFLOP/s, we
compute the theoretical floating-point operations for each kernel and matrix
and divide it by execution times of each tool.

The sequential implementation of SpMV CSR and SpTRSV CSR are used as a
baseline. Since the code for PIC is not publicly available, we created an in-
house inspector-executor implementation of their approach with feedback from
the authors of [2]. PIC was originally developed for a single thread, however,
we extended its support to be parallel and report the best performance between
the two implementations in all figures. A timeout of 4 hours was used for all
runs (including the inspector and the executor time). Thus, PIC does not have
data points in respective figures for large matrices. Also, our implementation
of ELL which is based on [13], has missing points in some figures because the
number of non-zero fill-ins exceeds the main memory of our test-bed
architecture.

###  V-A SpMV Performance

The performance of the SpMV kernel for DDF, MKL, CSR5, ELLPACK, baseline, and
PIC is shown in Figure 8. As shown, DDF is faster than all other tools for
over 88% of the matrices tested. DDF’s code is on average 1.93x, 1.42x,
10.57x, and 7.19x faster than in order MKL, CSR5, ELLPACK, and PIC
respectively on the Intel architecture. DDF’s code is on average 1.22x, 1.29x,
8.54x, and 10.93x faster than in order MKL, CSR5, ELLPACK, and PIC
respectively on the AMD architecture.

Figure 8: Top: Shows DDF speedups for the SpMV kernel over the single threaded
baseline SpMV CSR algorithm ran on an Intel architecture. Speedups also
reported for parallel MKL, parallel CSR5, parallel ELLPACK and single threaded
PIC [2] over the single threaded baseline SpMV CSR algorithm. Bottom: Shows
DDF, MKL, CSR5, ELLPACK and PIC speedups for the SpMV kernel using an AMD
architecture. Figure 9: Shows percentage of points for all matrices processed
by BLAS, PSC type 1 or PSC type 2 codelets for SpMV and SpTRSV.

To demonstrate the effect of partially strided codelets on the performance of
DDF, in Figure 9 we use a stacked bar for DDF. The stacks show the GFlop/s for
the baseline code (input to DDF), for DDF when it only mines for BLAS codelets
(DDF-baseline+BLAS and refer to DDF BLAS_only), and from running the entire
DDF algorithm, i.e. Algorithm 1, that mines for BLAS and PSC (shown with DDF
or DDF-baseline+BLAS+PSC in the figure). As shown, DDF is on average
10$\times$ faster than DDF BLAS_only which demonstrates the importance of
using PSC codelets. Figure 9 shows the percentage of operations that are
vectorized in the generated code from DDF with PSC I, PSC II, and BLAS
codelets, obtained by averaging over all matrices in the benchmark. As shown,
over 80% of the operations in SpMV are vectorized with PSC codelets.

The PSC codelets improve data locality in DDF’s generated code. Figure 10
shows the relation between DDF’s performance and the performance of the
parallel (OpenMP) version of the baseline code. Average memory access latency
[31] is used as a measure for locality and is computed by gathering the number
of misses and accesses to L1, L2, and LLC caches using the PAPI [32]
performance counters. Figure 10 shows the coefficient of determination or R2
is 0.83 which indicates a good correlation between speedup and the memory
access latency.

Figure 10: Shows the correlation between speedup over the baseline while
profiling the code with PAPI [32] and the relative memory cycle for matrices
used in the test set for SpMV. Relative memory cycle is the average memory
cycle of the baseline over the average memory cycle of the DDF framework
codes.

To further explain why DDF is faster than other tools, we conduct a few
experiments and report the resulting average over all matrices in this
paragraph. For MKL, we compute its average memory access latency, which is
4.36x slower than that of DDF, contributing to the worse performance compared
to DDF. To compare to CSR5, we count the number of instructions. CSR5 executes
1.73x more instructions compared to DDF. This is potentially due to the
overhead of the segmented sum calculations used in their approach to improving
vectorization and load balance. DDF is on average 10.57x faster than ELL,
because the percentage of fill-in in the ELL implementation is 800% relative
to the non-zero elements in the matrix, significantly increasing the number of
operations in their method. To conclude, the SpMV code of DDF is faster than
existing implementations because it improves data locality and/or reduces the
number of instructions via vectorization.

###  V-B SpTRSV Performance

Figure 11 compares the performance of SpTRSV using DDF, MKL, and Sympiler. The
stacked bar of DDF shows the additional performance gained from mining
partially strided codelets; the trend is similar to that in Figure 8 for SpMV.
On average DDF is faster than MKL, Sympiler, and the baseline 4.5$\times$,
1.79$\times$, and 3.62$\times$ respectively for the intel architecture. On the
AMD architecture, DDF is faster than MKL and Sympiler 2.66, 1.38$\times$
respectively.

As shown in Figure 11, partially strided codelets are the main contributors to
the overall performance of the DDF’s SpTRSV code. On average, DDF is
3.9$\times$ faster compared to when only BLAS codelets are generated. Similar
to SpMV, PSCs contribute to optimizing 78% of the operations over all matrices
(Figure 6b). However, the number of PSC II codelets has increased from 18% in
SpMV to 53% in SpTRSV. The number of computational regions with more than one
strided access function is small, due to the existing dependencies in the
SpTRSV kernel, thus, more PSC type II codelets are generated. Similar to the
SpMV kernel, the mined PSCs in DDF improve locality. The correlation
coefficient between the speedup and relative memory cycle is 0.67 which is
consistent with the trend in SpMV.

While MKL provides an efficient and vectorized implementation for single-
threaded SPTRV executions, it’s not optimized to execute on parallel
processors; the performance of MKL’s parallel code is similar to its serial
implementation. Sympiler performs well for matrices that contain row-blocks or
can be padded with up to 30% nonzeros to create row-blocks. These row-blocks
are converted to BLAS calls and thus improve locality. L-factors are generated
from a Cholesky factorization of the original matrices in the SuiteSparse
collection. Since L-factors are created with a supernodal factorization
process [33], the resulting matrices are well-structured and have many row-
blocks. As shown, Sympiler’s performance on L-factor matrices (IDs 31-50) is
close to that of DDF’s with an average speedup of 1.18$\times$, while for all
other matrices DDF’s performs on average time 2.1$\times$ better than
Sympiler.

Figure 11: Top: Shows breakdowns for DDF speedups, when enabling BLAS and PSC
codelet mining, MKL and Sympiler speedups over the single threaded baseline
SpTRSV CSR algorithm on Intel archtecture. Bottom: Shows DDF, MKL and Sympiler
speedups on AMD architecture.

###  V-C The Inspector Overhead

We compare the inspection time of DDF to that of other inspector-
executors/code-generation frameworks, i.e. Sympiler and PIC. We compute the
number of executor runs (NER) that amortize the cost of the inspector using
$\frac{\>\>Inspector\>\>Time}{Baseline\>Time-\>\>Executor\>\>Time}$. The
baseline time is obtained by running sequential implementation of the kernel.
The PIC inspector would timeout for over 83.3% of matrices because of its
large inspection overhead and code compilation time. For small matrices that
PIC would execute, more than one million executor runs are needed to amortize
the cost of the inspection. DDF’s inspection time is on average 0.5 seconds
with an average NER of 15 for SpMV. The inspection time of DDF and Sympiler
are similar with an average NER of less than 100 for both tools for SPTRV. The
largest matrices in our benchmark are inspected in less than 7 seconds with
DDF. Sparse kernels such as SpMV and SpTRSV are typically used in iterative
solvers, for example, to compute a residual in each iteration or to apply a
preconditioner per iteration. Even with preconditioning, these solvers
typically converge to a solution after tens of thousands of iterations [34,
35, 36] and hence inspector-executor frameworks such as DDF and Sympiler lead
to noticeable speedups as their inspection time overheads are amortized after
a few initial iterations of the solver.

##  VI Related Work

Numerous hand-optimized libraries [4, 37] and implementations [38, 39, 30, 40]
exist that optimize the performance of sparse matrix computations for
different parallel architectures and also optimize vectorization on a single
core. Libraries such as MKL and Eigen as well as implementations in [14, 41,
42, 43, 44] optimize the performance of SpMV on shared memory architectures
and improve SIMD vectorizability. A number of library implementations such
[13, 45, 46, 47, 11, 48, 12] reorganize data and computation to increase
opportunities for vectorization. A class of these libraries implement and
optimize sparse kernels based on available storage formats; for example [49]
optimizes SpMV based on ELLPACK. Other libraries work best for matrices
arising from specific applications such as [14], which optimizes SpMV for
large matrices from graph analytics, or KLU [50] which works best for circuit
simulation problems.

Domain-specific compilers use domain information to enable the application of
code transformations and optimizations such as vectorization. For example
domain-specific compilers such as [51, 52, 20] optimize signal processing
applications, stencil computations are optimized in [51, 52, 53], and matrix
assembly and mesh analysis is optimized in [54, 55]. While these compilers
provide highly-optimized code for the applications that they accelerate, they
do not generate specialized code for the input pattern and the computations do
not have indirect accesses.

Inspector-Executor approaches inspect the irregular access patterns of sparse
matrix computations at run-time to enable the automatic optimization of sparse
codes [56, 57, 58, 59, 23, 60, 61, 62, 63]. The index array accesses of the
sparse code is analyzed using an inspector and the information is used at run-
time to execute the code efficiently. The sparse polyhedral framework [64, 22,
65] uses uninterpreted function symbols to express regular and irregular
segments of sparse codes. As a result it is able to automatically generate
inspector executors at compile time that can resolve data dependencies in
sparse computations. These approaches do not generate code that is specialized
for the sparsity of the input matrix. Sympiler [1] and ParSy [28] are amongst
the inspector executor frameworks that inspect the matrix sparsity pattern and
as a result generate vectorized and parallel code specialized for the input
sparsity. Their optimizations for tiling and vactorization are based on
detecting row-blocks that primarily exist in matrices obtained from numerical
factorizations.

Augustine et al. [2] proposed an approach based on the Trace Reconstruction
Engine [66, 67] where polyhedral representations are built by inspecting the
sequence of addresses being accessed in the sparse matrix vector
multiplication. A followup to this work, proposes to use program guided
optimization for better vectorization [68]. These approaches lead to
generating code that is specialized for the sparsity pattern of the input
matrix and improves SIMD vectorization in SpMV. However, their work can only
support small matrices (below 0.5M nonzeros) because of inspector overheads
and because the size of the generated code increases with the matrix size.
Sparse computations with loop-carried dependencies such as the sparse
triangular solver are also not supported. DDF inspects memory address accesses
to find SIMD vectorizable code for sparse triangular solve and SpMV and
supports both small and large matrices.

##  VII Conclusion

In this work, we present partially strided codelets that enable the
vectorization of computation regions with unstrided memory accesses in sparse
matrix codes. We demonstrate how these codelets increase opportunities for
vectorization in sparse codes and also improve data locality in their
computation. A novel inspector-executor framework called DDF is proposed. DDF
uses an efficient inspector to mine for PSCs with a memory access
differentiation approach and as a result, generates highly efficient code for
sparse kernels. The performance of the DDF-generated code is compared to
state-of-the-art library implementations and other inspector-executor
frameworks for the sparse matrix-vector multiply and the sparse triangular
solve kernels. We also demonstrate that the inspection overhead of DDF is
negligible.

##  VIII Acknowledgments

This work was supported in part by NSERC Discovery Grants (RGPIN-06516,
DGECR00303), the Canada Research Chairs program, and U.S. NSF awards NSF
CCF-1814888, NSF CCF-1657175; used the Extreme Science and Engineering
Discovery Environment (XSEDE) [Towns et al. 2014] which is supported by NSF
grant number ACI-1548562; and was enabled in part by Compute Canada and
Scinet111www.computecanada.ca

## References

  * [1] K. Cheshmi, S. Kamil, M. M. Strout, and M. M. Dehnavi, “Sympiler: Transforming sparse matrix codes by decoupling symbolic analysis,” in _Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, SC 2017_. New York, NY, USA: Association for Computing Machinery, Inc, nov 2017, pp. 1–13. [Online]. Available: https://dl.acm.org/doi/10.1145/3126908.3126936
  * [2] T. Augustine, L. N. Pouchet, J. Sarma, and G. Rodríguez, “Generating piecewise-regular code from irregular structures,” in _Proceedings of the ACM SIGPLAN Conference on Programming Language Design and Implementation (PLDI)_. New York, NY, USA: Association for Computing Machinery, jun 2019, pp. 625–639. [Online]. Available: https://dl.acm.org/doi/10.1145/3314221.3314615
  * [3] A. Buluç, J. T. Fineman, M. Frigo, J. R. Gilbert, and C. E. Leiserson, “Parallel sparse matrix-vector and matrix-transpose-vector multiplication using compressed sparse blocks,” in _Proceedings of the twenty-first annual symposium on Parallelism in algorithms and architectures_ , 2009, pp. 233–244.
  * [4] _Intel Math Kernel Library. Reference Manual_. Intel Corporation, 2009, santa Clara, USA. ISBN 630813-054US.
  * [5] M. Kreutzer, G. Hager, G. Wellein, H. Fehske, and A. R. Bishop, “A unified sparse matrix data format for efficient general sparse matrix-vector multiply on modern processors with wide SIMD units,” _SIAM Journal on Scientific Computing_ , vol. 36, no. 5, pp. C401–C423, jul 2013. [Online]. Available: http://arxiv.org/abs/1307.6209http://dx.doi.org/10.1137/130930352
  * [6] X. Liu, M. Smelyanskiy, E. Chow, and P. Dubey, “Efficient sparse matrix-vector multiplication on x86-based many-core processors,” in _Proceedings of the International Conference on Supercomputing_. New York, New York, USA: ACM Press, 2013, pp. 273–282. [Online]. Available: http://dl.acm.org/citation.cfm?doid=2464996.2465013
  * [7] W. Liu and B. Vinter, “CSR5: An Efficient Storage Format for Cross-Platform Sparse Matrix-Vector Multiplication,” _Proceedings of the International Conference on Supercomputing_ , vol. 2015-June, pp. 339–350, mar 2015. [Online]. Available: http://arxiv.org/abs/1503.05032
  * [8] W. T. Tang, R. Zhao, M. Lu, Y. Liang, H. P. Huyng, X. Li, and R. S. M. Goh, “Optimizing and auto-tuning scale-free sparse matrix-vector multiplication on Intel Xeon Phi,” in _Proceedings of the 2015 IEEE/ACM International Symposium on Code Generation and Optimization, CGO 2015_. Institute of Electrical and Electronics Engineers Inc., mar 2015, pp. 136–145.
  * [9] L. Chen, P. Jiang, and G. Agrawal, “Exploiting recent SIMD architectural advances for irregular applications,” in _Proceedings of the 14th International Symposium on Code Generation and Optimization, CGO 2016_. New York, NY, USA: Association for Computing Machinery, Inc, feb 2016, pp. 47–58. [Online]. Available: https://dl.acm.org/doi/10.1145/2854038.2854046
  * [10] B. Xie, W. Gao, J. Zhan, Z. Jia, L. Zhang, X. Liu, and X. He, “CVR: Efficient Vectorization of SpMV on X86 Processors,” in _CGO 2018 - Proceedings of the 2018 International Symposium on Code Generation and Optimization_ , vol. 2018-February. New York, NY, USA: Association for Computing Machinery, Inc, feb 2018, pp. 149–162. [Online]. Available: https://dl.acm.org/doi/10.1145/3168818
  * [11] W. Liu and B. Vinter, “Csr5: An efficient storage format for cross-platform sparse matrix-vector multiplication,” in _Proceedings of the 29th ACM on International Conference on Supercomputing_ , 2015, pp. 339–350.
  * [12] B. Xie, J. Zhan, X. Liu, W. Gao, Z. Jia, X. He, and L. Zhang, “Cvr: Efficient vectorization of spmv on x86 processors,” in _Proceedings of the 2018 International Symposium on Code Generation and Optimization_ , 2018, pp. 149–162.
  * [13] Y. Saad, _Iterative methods for sparse linear systems_. SIAM, 2003.
  * [14] S. Yesil, A. Heidarshenas, A. Morrison, and J. Torrellas, “Speeding up spmv for power-law graph analytics by enhancing locality vectorization,” in _International Conference for High Performance Computing, Networking, Storage and Analysis, SC_ , vol. 2020-November. IEEE Computer Society, nov 2020, pp. 1–15.
  * [15] F. Vazquez, G. Ortega, J.-J. Fernández, and E. M. Garzón, “Improving the performance of the sparse matrix vector product with gpus,” in _2010 10th IEEE International Conference on Computer and Information Technology_. IEEE, 2010, pp. 1146–1151.
  * [16] P. Boulet and P. Feautrier, “Scanning polyhedra without Do-loops,” in _Parallel Architectures and Compilation Techniques - Conference Proceedings, PACT_. Institute of Electrical and Electronics Engineers Inc., 1998, pp. 4–11.
  * [17] C. Chen, “Polyhedra scanning revisited,” in _Proceedings of the ACM SIGPLAN Conference on Programming Language Design and Implementation (PLDI)_. New York, New York, USA: ACM Press, 2012, pp. 499–508. [Online]. Available: http://dl.acm.org/citation.cfm?doid=2254064.2254123
  * [18] D. S. Kershaw, “The incomplete Cholesky-conjugate gradient method for the iterative solution of systems of linear equations,” pp. 43–65, jan 1978.
  * [19] M. Püschel, J. M. Moura, J. R. Johnson, D. Padua, M. M. Veloso, B. W. Singer, J. Xiong, F. Franchetti, A. Gačić, Y. Voronenko, K. Chen, R. W. Johnson, and N. Rizzolo, “SPIRAL: Code generation for DSP transforms,” in _Proceedings of the IEEE_ , vol. 93, no. 2. Institute of Electrical and Electronics Engineers Inc., 2005, pp. 232–273. [Online]. Available: https://experts.illinois.edu/en/publications/spiral-code-generation-for-dsp-transforms
  * [20] D. G. Spampinato and M. Püschel, “A basic linear algebra compiler,” in _Proceedings of the 12th ACM/IEEE International Symposium on Code Generation and Optimization, CGO 2014_. New York, NY, USA: Association for Computing Machinery, feb 2014, pp. 23–32. [Online]. Available: http://dl.acm.org/doi/10.1145/2544137.2544155
  * [21] A. Tiwari, C. Chen, J. Chame, M. Hall, and J. K. Hollingsworth, “A scalable auto-tuning framework for compiler optimization,” in _IPDPS 2009 - Proceedings of the 2009 IEEE International Parallel and Distributed Processing Symposium_ , 2009.
  * [22] M. M. Strout, G. Georg, and C. Olschanowsky, “Set and relation manipulation for the sparse polyhedral framework,” in _International Workshop on Languages and Compilers for Parallel Computing_. Springer, 2012, pp. 61–75.
  * [23] A. Venkat, M. Hall, and M. Strout, “Loop and data transformations for sparse matrix code,” in _Proceedings of the ACM SIGPLAN Conference on Programming Language Design and Implementation (PLDI)_ , vol. 2015-June. New York, NY, USA: Association for Computing Machinery, jun 2015, pp. 521–532. [Online]. Available: https://dl.acm.org/doi/10.1145/2737924.2738003
  * [24] M. M. Strout, M. Hall, and C. Olschanowsky, “The sparse polyhedral framework: Composing compiler-generated inspector-executor code,” _Proceedings of the IEEE_ , vol. 106, no. 11, pp. 1921–1934, 2018.
  * [25] L. S. Blackford, A. Petitet, R. Pozo, K. Remington, R. C. Whaley, J. Demmel, J. Dongarra, I. Duff, S. Hammarling, G. Henry _et al._ , “An updated set of basic linear algebra subprograms (blas),” _ACM Transactions on Mathematical Software_ , vol. 28, no. 2, pp. 135–151, 2002.
  * [26] T. A. Davis and Y. Hu, “The university of florida sparse matrix collection,” _ACM Transactions on Mathematical Software (TOMS)_ , vol. 38, no. 1, pp. 1–25, 2011.
  * [27] ——, “The University of Florida Sparse Matrix Collection,” _ACM Transactions on Mathematical Software_ , vol. 38, no. 1, pp. 1–25, nov 2011. [Online]. Available: https://dl.acm.org/doi/10.1145/2049662.2049663
  * [28] K. Cheshmi, S. Kamil, M. M. Strout, and M. M. Dehnavi, “ParSy: Inspection and transformation of sparse matrix computations for parallelism,” in _Proceedings - International Conference for High Performance Computing, Networking, Storage, and Analysis, SC 2018_. Institute of Electrical and Electronics Engineers Inc., mar 2019, pp. 779–793.
  * [29] E. Wang, Q. Zhang, B. Shen, G. Zhang, X. Lu, Q. Wu, and Y. Wang, _Intel Math Kernel Library_ , 05 2014, pp. 167–188.
  * [30] T. A. Davis and W. Hager, “Cholmod: supernodal sparse cholesky factorization and update/downdate,” 2005.
  * [31] J. Hennessy and D. Patterson, _Computer Architecture - A Quantitative Approach_ , 01 2007.
  * [32] D. Terpstra, H. Jagode, H. You, and J. Dongarra, “Collecting performance data with papi-c,” in _Tools for High Performance Computing 2009_. Springer, 2010, pp. 157–173.
  * [33] Y. Chen, T. A. Davis, W. W. Hager, and S. Rajamanickam, “Algorithm 887: Cholmod, supernodal sparse cholesky factorization and update/downdate,” _ACM Transactions on Mathematical Software (TOMS)_ , vol. 35, no. 3, pp. 1–14, 2008.
  * [34] M. Benzi, J. K. Cullum, and M. Tuma, “Robust approximate inverse preconditioning for the conjugate gradient method,” _SIAM Journal on Scientific Computing_ , vol. 22, no. 4, pp. 1318–1332, 2000.
  * [35] D. S. Kershaw, “The incomplete cholesky-conjugate gradient method for the iterative solution of systems of linear equations,” _Journal of computational physics_ , vol. 26, no. 1, pp. 43–65, 1978.
  * [36] M. Papadrakakis and N. Bitoulas, “Accuracy and effectiveness of preconditioned conjugate gradient algorithms for large and ill-conditioned problems,” _Computer methods in applied mechanics and engineering_ , vol. 109, no. 3-4, pp. 219–232, 1993.
  * [37] G. Guennebaud, B. Jacob _et al._ , “Eigen v3,” http://eigen.tuxfamily.org, 2010.
  * [38] X. S. Li and J. W. Demmel, “SuperLU_DIST: A scalable distributed-memory sparse direct solver for unsymmetric linear systems,” _ACM Transactions on Mathematical Software_ , vol. 29, no. 2, pp. 110–140, jun 2003\. [Online]. Available: https://dl.acm.org/doi/10.1145/779359.779361
  * [39] S. A. Hutchinson, J. N. Shadid, and R. S. Tuminaro, “Aztec User’s Guide Version 1.1,” Tech. Rep., 1995.
  * [40] T. A. Davis, “Algorithm 832: Umfpack v4. 3—an unsymmetric-pattern multifrontal method,” _ACM Transactions on Mathematical Software (TOMS)_ , vol. 30, no. 2, pp. 196–199, 2004.
  * [41] R. Li and Y. Saad, “Gpu-accelerated preconditioned iterative linear solvers,” _The Journal of Supercomputing_ , vol. 63, no. 2, pp. 443–466, 2013.
  * [42] S. Williams, L. Oliker, R. Vuduc, J. Shalf, K. Yelick, and J. Demmel, “Optimization of sparse matrix-vector multiplication on emerging multicore platforms,” in _SC’07: Proceedings of the 2007 ACM/IEEE Conference on Supercomputing_. IEEE, 2007, pp. 1–12.
  * [43] S. Kamin, M. J. Garzarán, B. Aktemur, D. Xu, B. Yılmaz, and Z. Chen, “Optimization by runtime specialization for sparse matrix-vector multiplication,” in _Proceedings of the 2014 International Conference on Generative Programming: Concepts and Experiences_ , 2014, pp. 93–102.
  * [44] D. Merrill and M. Garland, “Merge-based parallel sparse matrix-vector multiplication,” in _SC’16: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_. IEEE, 2016, pp. 678–689.
  * [45] R. W. Vuduc, _Automatic performance tuning of sparse matrix kernels_. University of California, Berkeley, 2003.
  * [46] R. W. Vuduc and H.-J. Moon, “Fast sparse matrix-vector multiplication by exploiting variable block structure,” in _International Conference on High Performance Computing and Communications_. Springer, 2005, pp. 807–816.
  * [47] Y. Li, P. Xie, X. Chen, J. Liu, B. Yang, S. Li, C. Gong, X. Gan, and H. Xu, “Vbsf: a new storage format for simd sparse matrix–vector multiplication on modern processors,” _The Journal of Supercomputing_ , vol. 76, no. 3, pp. 2063–2081, 2020.
  * [48] W. T. Tang, R. Zhao, M. Lu, Y. Liang, H. P. Huyng, X. Li, and R. S. M. Goh, “Optimizing and auto-tuning scale-free sparse matrix-vector multiplication on intel xeon phi,” in _2015 IEEE/ACM International Symposium on Code Generation and Optimization (CGO)_. IEEE, 2015, pp. 136–145.
  * [49] M. Kreutzer, G. Hager, G. Wellein, H. Fehske, and A. R. Bishop, “A unified sparse matrix data format for efficient general sparse matrix-vector multiplication on modern processors with wide simd units,” _SIAM Journal on Scientific Computing_ , vol. 36, no. 5, pp. C401–C423, 2014.
  * [50] T. A. Davis and E. Palamadai Natarajan, “Algorithm 907: Klu, a direct sparse solver for circuit simulation problems,” _ACM Transactions on Mathematical Software (TOMS)_ , vol. 37, no. 3, pp. 1–17, 2010.
  * [51] N. D. Jones, C. K. Gomard, and P. Sestoft, _Partial evaluation and automatic program generation_. Peter Sestoft, 1993.
  * [52] F. Quilleré, S. Rajopadhye, and D. Wilde, “Generation of efficient nested loops from polyhedra,” _International journal of parallel programming_ , vol. 28, no. 5, pp. 469–498, 2000.
  * [53] D. G. Spampinato and M. Püschel, “A basic linear algebra compiler,” in _Proceedings of Annual IEEE/ACM International Symposium on Code Generation and Optimization_ , 2014, pp. 23–32.
  * [54] M. S. Alnæs, A. Logg, K. B. Ølgaard, M. E. Rognes, and G. N. Wells, “Unified form language: A domain-specific language for weak formulations of partial differential equations,” _ACM Transactions on Mathematical Software (TOMS)_ , vol. 40, no. 2, pp. 1–37, 2014.
  * [55] F. Luporini, D. A. Ham, and P. H. Kelly, “An algorithm for the optimization of finite element integration loops,” _ACM Transactions on Mathematical Software (TOMS)_ , vol. 44, no. 1, pp. 1–26, 2017.
  * [56] J. Shin, M. W. Hall, J. Chame, C. Chen, P. F. Fischer, and P. D. Hovland, “Speeding up Nek5000 with autotuning and specialization,” in _Proceedings of the International Conference on Supercomputing_. New York, New York, USA: ACM Press, 2010, pp. 253–262. [Online]. Available: http://portal.acm.org/citation.cfm?doid=1810085.1810120
  * [57] Y. Tang, R. A. Chowdhury, B. C. Kuszmaul, C. K. Luk, and C. E. Leiserson, “The pochoir stencil compiler,” in _Annual ACM Symposium on Parallelism in Algorithms and Architectures_. New York, New York, USA: ACM Press, 2011, pp. 117–128. [Online]. Available: http://portal.acm.org/citation.cfm?doid=1989493.1989508
  * [58] H. L. A. van der Spek and H. A. G. Wijshoff, “Sublimation: Expanding data structures to enable data instance specific optimizations,” in _Languages and Compilers for Parallel Computing_ , K. Cooper, J. Mellor-Crummey, and V. Sarkar, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2011, pp. 106–120.
  * [59] N. Vasilache, C. Bastoul, and A. Cohen, “Polyhedral code generation in the real world,” in _International Conference on Compiler Construction_. Springer, 2006, pp. 185–201.
  * [60] G. Agrawal, J. Saltz, and R. Das, “Interprocedural partial redundancy elimination and its application to distributed memory compilation,” _ACM SIGPLAN Notices_ , vol. 30, no. 6, pp. 258–269, 1995.
  * [61] R. Das, P. Havlak, J. Saltz, and K. Kennedy, “Index array flattening through program transformation,” in _Proceedings of the 1995 ACM/IEEE conference on Supercomputing_ , 1995, pp. 70–es.
  * [62] G. Rodríguez and L. Pouchet, “Polyhedral Modeling of Immutable Sparse Matrices,” _impact.gforge.inria.fr_. [Online]. Available: http://impact.gforge.inria.fr/impact2018/papers/modeling-immutable-sparsemat.pdf
  * [63] J. Saltz, K. Crowley, R. Michandaney, and H. Berryman, “Run-time scheduling and execution of loops on message passing machines,” _Journal of Parallel and Distributed Computing_ , vol. 8, no. 4, pp. 303–312, 1990.
  * [64] A. LaMielle and M. M. Strout, “Enabling code generation within the sparse polyhedral framework,” _Technical report, Technical Report CS-10-102_ , 2010\.
  * [65] A. Venkat, M. S. Mohammadi, J. Park, H. Rong, R. Barik, M. M. Strout, and M. Hall, “Automating wavefront parallelization for sparse matrix computations,” in _SC’16: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_. IEEE, 2016, pp. 480–491.
  * [66] G. Rodríguez, J. M. Andión, M. T. Kandemir, and J. Touriño, “Trace-based affine reconstruction of codes,” in _Proceedings of the 14th International Symposium on Code Generation and Optimization, CGO 2016_. New York, NY, USA: Association for Computing Machinery, Inc, feb 2016, pp. 139–149. [Online]. Available: https://dl.acm.org/doi/10.1145/2854038.2854056
  * [67] A. Ketterlin and P. Clauss, “Prediction and trace compression of data access addresses through nested loop recognition,” in _Proceedings of the 2008 CGO - Sixth International Symposium on Code Generation and Optimization_. New York, New York, USA: ACM Press, 2008, pp. 94–103. [Online]. Available: http://portal.acm.org/citation.cfm?doid=1356058.1356071
  * [68] M. Selva, F. Gruber, D. Sampaio, C. Guillon, L. N. Pouchet, and F. Rastello, “Building a polyhedral representation from an instrumented execution: Making dynamic analyses of nonaffine programs scalable,” _ACM Transactions on Architecture and Code Optimization_ , vol. 16, no. 4, pp. 1–26, dec 2019. [Online]. Available: https://dl.acm.org/doi/10.1145/3363785