11institutetext: Technical University Munich, Germany 22institutetext:
Helmholtz Munich, Germany 33institutetext: Munich Center of Machine Learning
(MCML), Germany 44institutetext: Universitat de Barcelona, Spain
55institutetext: King’s College London, UK

# Progressive Growing of Patch Size: Resource-Efficient Curriculum Learning
for Dense Prediction Tasks

Stefan M. Fischer 112233 Lina Felsner 1122 Richard Osuala 112244 Johannes
Kiechle 112233 Daniel M. Lang 1122 Jan C. Peeken 1122 Julia A. Schnabel
11223355

###### Abstract

In this work, we introduce Progressive Growing of Patch Size, a resource-
efficient implicit curriculum learning approach for dense prediction tasks.
Our curriculum approach is defined by growing the patch size during model
training, which gradually increases the task’s difficulty. We integrated our
curriculum into the nnU-Net framework and evaluated the methodology on all 10
tasks of the Medical Segmentation Decathlon. With our approach, we are able to
substantially reduce runtime, computational costs, and CO2 emissions of
network training compared to classical constant patch size training. In our
experiments, the curriculum approach resulted in improved convergence. We are
able to outperform standard nnU-Net training, which is trained with constant
patch size, in terms of Dice Score on 7 out of 10 MSD tasks while only
spending roughly 50% of the original training runtime. To the best of our
knowledge, our Progressive Growing of Patch Size is the first successful
employment of a sample-length curriculum in the form of patch size in the
field of computer vision. Our code is publicly available at
https://github.com/compai-lab/2024-miccai-fischer.

###### Keywords:

Segmentation nnU-Net Curriculum Learning Resource Efficiency Medical
Segmentation Decathlon

##  1 Introduction

Automatic medical image segmentation is dominated by deep learning based
methods. Most research focuses on the development of new architectural
concepts, introducing convolution-based [15, 7], transformer-based [3, 13] or
hybrid approaches [5, 4] to improve downstream segmentation performance. In
contrast, few works in the medical imaging domain focus on the actual training
process of networks, and models are still trained mainly by random training
data sampling. Inspired by humans, Bengio et al. [2] introduced the concept of
curriculum learning, teaching models to first solve simple tasks and,
subsequent to that, concentrate on harder tasks. They showed that ordering
samples from easy-to-hard speeds up convergence and can, therefore, also
improve performance compared to training with random sampling. Those training
techniques can reduce the high costs of model training in terms of time,
energy, and carbon footprint [17].

Different approaches have been used in the image processing field to establish
such a curriculum sample ordering. Human annotations or expert knowledge can
be included as problem-specific measures [2, 10, 19]. Task difficulty can also
be directly linked to sample class membership, which was utilized for fracture
classification [10]. Another frequently applied measure in the medical domain
is the inter-rater expert agreement, which considerably increases annotation
costs [10, 19].

A more universal approach is rating the task difficulty automatically, or by
incorporating processes that increase the task difficulty synthetically.
Sample difficulty estimation can directly be computed from the network’s
training sample loss. Jesson et al. [9] utilized hard-negative mining and an
adaptive oversampling scheme for lung cancer segmentation. Missing modalities
in MRI processing have been the motivation for Havei et al. [6] to train a
brain tumor segmentation model that is robust against such missing data
scenarios. In their curriculum learning approach, they randomly drop MRI
channels with increasing probability. In the field of natural image
generation, Kerras et al. [11] defined an implicit curriculum by progressively
growing their image generator and discriminator layer-wise during training.
With each additional layer, the output resolution is doubled, resulting in
more difficult tasks. This increased convergence speed and reduced training
runtime while achieving higher model performance compared to standard
training. Zhao et al. [22] directly transferred this idea to semantic
segmentation for their task of cervical nuclei segmentation using a U-net
architecture, which grows starting from the bottleneck.

Task difficulty can also be directly defined by sample-length, which is
already used in the natural language processing field in the form of sentence
length [18, 21, 14]. In contrast, in the computer vision domain, the sample
length, which would refer to image size or patch size, has, to the best of our
knowledge, not used yet. An indirect utilization of image size is used in the
Progressive Growing of GANs and its adaption to segmentation, but changes in
image size are a result of the addition of network layers [11, 22]. In
contrast to that growing image resolution curriculum that would start
segmenting a lesion from a low-resolution version of the full image, we assume
that starting with segmenting a foreground object within a smaller region-of-
interest patch is easier. This was our motivation to build our curriculum on
image size instead of image resolution. In this work, we successfully apply
the curriculum approach of growing patch size. Our main contributions in this
work are:

  1. 1.

Introduction of Progressive Growing of Patch Size curriculum for semantic
segmentation and integration into the nnU-Net framework

  2. 2.

Empirical verification of reduced computational costs and improved convergence
compared to classical constant patch size training

  3. 3.

Validation of the robustness of our approach on the 10 tasks of the Medical
Segmentation Decathlon (MSD) [1]

##  2 Methods and Materials

###  2.1 Progressive Growing of Patch Size as Curriculum

We establish our curriculum by changing the patch size of the input volume. We
assume training on smaller patches is an easier task than full images or
volumes, as the frequency of foreground pixels decreases in larger patches.
Furthermore, growing the patch size also results in an increased global
context that can be inferred from that patch.

For our approach, we start training the network with a minimal patch size and
then linearly increase the patch size to a maximal patch size during training,
training at each patch size for the same number of iterations. During the
model inference, the maximal possible patch size is used. A sketch of our
curriculum is shown in Fig. 1.

Figure 1: Proposed Progressive Growing of Patch Size curriculum. Our
curriculum is illustrated for the lung cancer segmentation with cancer
bounding boxes in yellow. A fully convolutional network is able to handle
inputs of different sizes. Training with our curriculum starts by training the
network with minimal patch size and, with training progress, increasing the
patch size until the final patch size is reached. The ratio between foreground
and background voxels is bigger for small patch sizes and decreases for bigger
patch sizes. In contrast, the global context that can be inferred from the
patch is growing with the patch size. For inference, the maximal patch size is
used.

###  2.2 Implementation with nnU-Net

Fully convolutional networks are able to process inputs of varying sizes,
which is the basis for our approach. While the GPU memory restricts the
maximal patch size we can set, the network architecture itself restricts the
minimal patch sizes that can be processed. Starting from the smallest possible
patch size, we increase the patch size in the smallest possible steps to keep
the similarity between different patch size stages as close as possible.

For our experiments, we use the current state-of-the-art medical segmentation
network nnU-Net, which follows the U-net architecture [15] and is a fully
self-configuring pipeline, as detailed in [7]. The experiments were performed
on the 3D patch-based segmentation version of nnU-Net (3d_fullres), which has
shown, on average, the best performances in the MSD in [8]. Models trained
with the standard training of the nnU-Net are referred to as models trained
with the Constant Patch Size Training Scheme (CPS). For our Progressive
Growing of Patch Size Training Scheme (PGPS), we only change the patch size
compared to CPS.

###  2.3 Utilizing Smaller Patch Sizes for Bigger Batch Sizes

Training a network with a smaller patch size reduces GPU memory consumption.
Thus, PGPS has lower GPU memory consumption during most of the training
compared to CPS. This enables the use of larger batch sizes in the lower patch
size phases of PGPS at low GPU memory costs. With a higher batch size, we
expect an increase in performance and convergence speed. Besides, we force a
foreground-background ratio of one, which is the default effective ratio of
nnU-Net that uses a batch size of two, a result of it’s maximal patch size
heuristic.

##  3 Experiments and Results

###  3.1 Lung Cancer Segmentation

####  3.1.1 Experiments

We compare nnU-Net instances trained via CPS (default nnU-Net) and nnU-Net
instances trained via our PGPS on the MSD task of lung cancer segmentation. To
evaluate the effect of growing patch size, we also train models with Random
Patch Size Sampling (RPSS). For that, at each training iteration, a random
patch size of all patch sizes used for PGPS is picked, and a training
iteration with the chosen patch size is performed. At testing time, all
configurations use the same patch size, which is the nnU-Net default patch
size.

To evaluate the convergence properties of the three different curricula, we
train models for different numbers of training iterations per epoch. The
standard nnU-Net 250 iterations per epoch refers to 100% scenario. Moreover,
we also train models with only 10%, 25%, 50% of iterations per epoch. For each
model training, we track the average training runtime and compute the mean
Dice Score and the general number of voxels shown to the network during
training.

Performance regarding training iterations, iterated voxels, and overall
training runtime for trained models are plotted in Fig. 2. Exact values are
given in the Supplemental Table 3. Performance of the 100% iterations per
epoch training of the CPS are taken from [8].

Figure 2: Segmentation performance of models trained with Progressive Growing of Patch Size (PGPS), Random Patch Size Sampling (RPSS) and Constant Patch Size (CPS) for different numbers of training iterations per epoch on the MSD Lung Cancer Task. Models were trained with 10%, 25%, 50%, and 100% of the training iterations per epoch, while 100% represents the default 250 iterations per epoch of nnU-Net. Dice Scores are averaged over a 5-fold cross-validation. All three plots are referring to the same models but are evaluated regarding different measures. Table 1: Performance of Constant Patch Size Training Scheme (CPS), Progressive Growing of Patch Size Training Scheme (PGPS), and PGPS with increased batch size (PGPS+) on Medical Segmentation Decathlon. CPS refers to standard nnU-Net training. (Dice Score: Evaluated in 5-fold Cross-Validation as in [8]; P-value: One-sided paired T-test on validation Dice Scores of all samples against CPS (* refers to significant difference [P-value $<$ 0.05]); Runtime: Training runtime per fold on one NVIDIA A100-SXM4-80GB GPU; Voxels Shown: total number of voxels iterated during training normalized to CPS training; CO2-eq: mean CO2-equivalent for training a single nnU-Net instance) MSD Task |  Scheme |  Dice Score |  P-value |  Runtime |  Voxels Shown |  CO2-eq
---|---|---|---|---|---|---
|  CPS |  0.7411 |  -* |  27.53 h |  100.00 % |  10.44 kg
Brain |  PGPS |  0.7412 |  0.5972* |  11.29 h |  38.08 % |  4.31 kg
|  PGPS+ |  0.7421 |  0.1173* |  13.83 h |  51.78 % |  5.21 kg
|  CPS |  0.9328 |  -* |  15.04 h |  100.00 % |  6.24 kg
Heart |  PGPS |  0.9321 |  0.2567* |  6.44 h |  35.35 % |  2.64 kg
|  PGPS+ |  0.9328 |  0.4763* |  7.06 h |  39.35 % |  2.89 kg
|  CPS |  0.7971 |  -* |  11.09 h |  100.00 % |  4.70 kg
Liver |  PGPS |  0.7891 |  0.0992* |  5.14 h |  38.08 % |  2.14 kg
|  PGPS+ |  0.7938 |  0.2369* |  6.70 h |  51.78 % |  2.77 kg
|  CPS |  0.8891 |  -* |  2.14 h |  100.00 % |  0.90 kg
Hippocampus |  PGPS |  0.8911 |  0.0073* |  1.37 h |  35.05 % |  0.54 kg
|  PGPS+ |  0.8907 |  0.0343* |  1.42 h |  35.23 % |  0.56 kg
|  CPS |  0.7537 |  -* |  10.14 h |  100.00 % |  4.33 kg
Prostate |  PGPS |  0.7566 |  0.4021* |  5.18 h |  30.32 % |  2.11 kg
|  PGPS+ |  0.7531 |  0.4679* |  5.75 h |  36.33 % |  2.31 kg
|  CPS |  0.7211 |  -* |  13.55 h |  100.00 % |  5.59 kg
Lung |  PGPS |  0.7263 |  0.3163* |  5.70 h |  35.35 % |  2.35 kg
|  PGPS+ |  0.7333 |  0.1484* |  6.13 h |  39.35 % |  2.55 kg
|  CPS |  0.6745 |  -* |  11.01 h |  100.00 % |  4.68 kg
Pancreas |  PGPS |  0.6824 |  0.0681* |  5.11 h |  33.93 % |  2.07 kg
|  PGPS+ |  0.6822 |  0.0738* |  5.19 h |  35.92 % |  2.30 kg
|  CPS |  0.6837 |  -* |  15.12 h |  100.00 % |  6.06 kg
Hepatic Vessel |  PGPS |  0.6782 |  0.0210* |  6.18 h |  36.19 % |  2.51 kg
|  PGPS+ |  0.6871 |  0.0633* |  6.65 h |  40.36 % |  2.72 kg
|  CPS |  0.9638 |  -* |  11.24 h |  100.00 % |  4.76 kg
Spleen |  PGPS |  0.9621 |  0.3803* |  5.22 h |  36.30 % |  2.15 kg
|  PGPS+ |  0.9654 |  0.4147* |  5.37 h |  45.81 % |  2.31 kg
|  CPS |  0.4553 |  -* |  10.94 h |  100.00 % |  3.77 kg
Colon |  PGPS |  0.4925 |  0.0046* |  4.36 h |  32.21 % |  1.69 kg
|  PGPS+ |  0.4967 |  0.0087* |  4.50 h |  34.55 % |  1.77 kg

####  3.1.2 Results

Focusing on the number of training iterations per epoch (Fig. 2 [left]), PGPS
models outperformed RPSS and CPS models in Dice Score for each trained
configuration. RPSS is outperformed by CPS and PGPS for all training
configurations except for 50%, where it outperforms CPS. When comparing the
same models regarding training efficiency, in terms of shown voxels during
training (Fig. 2 [center]) and training runtime (Fig. 2 [right]), PGPS and
RPSS have drastically increased efficiency. The number of shown training
voxels for PGPS and RPSS is substantially reduced for each configuration to
only 35.5% of the original CPS training amount. Training runtime is reflecting
the trend of voxels shown. Training with 100% iterations per epoch is only
running for 5.7 hours for PGPS instead of 13.55 hours for CPS.

###  3.2 Medical Segmentation Decathlon

####  3.2.1 Experiments

To evaluate the robustness of PGPS on a variety of medical segmentation tasks,
we tested our curriculum learning approach on all ten different tasks of the
MSD. Detailed descriptions of the tasks are given in [1]. Various anatomies
have to be segmented, including tumors, vessels, and healthy organs on single
or multi-modality input, covering CT and MRI images. The MSD covers binary and
multi-class segmentation tasks and dataset sample sizes ranging from 30 to 750
samples. Dice Score, runtime, and total number of voxels shown during training
are computed over a 5-fold cross-validation as in [8]. Besides, we adopt
CodeCarbon [16] to track the CO2-equivalents of our training process. Dice
Scores of CPS models are taken from [8]. Because of that, we extrapolated the
runtime values for CPS from the runtimes of PGPS by multiplying the averaged
measured time for one epoch of the maximal patch size times the number of
training epochs. Used nnU-Net model architectures and training hyperparameters
can be seen for each MSD task in the Appendix of [8]. Concrete patch sizes
used for PGPS are given in the Supplemental Table 3. Furthermore, we also
compute P-values for paired one-sided T-tests to assess if the performance
differs significantly. For that, we pair validation Dice Scores of each sample
over all cross-validation splits regarding CPS and PGPS.

Besides training with PGPS, we also repeat experiments with increased batch
size. In our approach, we increase the batch size of the current patch size
stage to the extent that the total number of voxels in an input tensor is
equal to or lower than that of the following patch size stage. In this way, we
avoid drops in the number of input voxels per iteration. We refer to this
adaption of PGPS as PGPS+ in the later text. The concrete patch and batch
sizes used for PGPS and PGPS+ are given in the Supplemental Table 3.

####  3.2.2 Results

Results on the MSD tasks are given in Table 1. In Fig. 3, we plotted the fold-
wise Dice Score differences between our new training curricula and standard
CPS. On averaged Dice Scores over the 5-fold cross-validation, the PGPS models
outperform CPS models in 6 out of 10 tasks. In the hippocampus and the colon
cancer task, PGPS outperforms CPS significantly on the sample-based T-Test
(P-value $<$ 0.05), while PGPS is outperformed significantly by CPS on the
hepatic vessel task. When PGPS outperformed CPS, it achieved an average fold-
wise increase of $1.8\%\pm 2.9\%$ in Dice Score compared to CPS, while for
tasks CPS outperformed PGPS, there was only an average fold-wise difference of
$0.5\%\pm 0.4\%$ in performance. Full training of the nnU-Net with our PGPS
curriculum reduces the runtime on average to $46.09\%\pm 6.8\%$ compared to
CPS (default nnU-Net training) by only iterating $35.09\%\pm 2.30\%$ of the
CPS’ seen voxels and reducing CO2 equivalents to $48.22\%\pm 9.34\%$.

Figure 3: Performance of Progressive Growing of Patch Size (PGPS), PGPS with
increased batch size (PGPS+) and Constant Patch Size (CPS) on the MSD tasks.
Here, we show the fold-wise mean difference in Dice Score and its standard
error of the 5-fold cross-validation. A positive difference represents wins of
the PGPS curricula over the CPS training (here marked in dark blue). On this
measure, PGPS outperforms CPS in 6 of 10 tasks, and PGPS+ outperforms CPS in 8
of 10 tasks.

With the additional adjustment of batch size for PGPS+, we are able to
outperform CPS in 7 out of 10 MSD tasks on the mean of fold-wise Dice Scores.
Sample-based significant performance increases were observed for the
hippocampus, pancreas, and colon cancer tasks. When PGPS+ outperformed CPS, it
achieved an average increase of $1.7\%\pm 3.1\%$ in fold-wise average Dice
Score compared to CPS, while when it lost, there was only an average
difference of $0.2\%\pm 0.2\%$. Comparing the performance differences between
each fold, our PGPS+ even outperforms CPS in 8 out of 10 tasks, which can be
seen in Fig. 3. The model training only takes up an average of $50.59\%\pm
7.6\%$ of the original time of CPS, while iterated voxels are reduced to
$41.05\%\pm 6.20\%$ and CO2 equivalents to $50.59\%\pm 5.55\%$.

##  4 Discussion and Conclusion

We introduced a novel curriculum learning approach based on increasing the
patch size during training and compared it with constant patch size training.

Our curriculum approach showed improved convergence in the experiments, as it
significantly outperforms classical CPS on the MSD most difficult task of
colon cancer by a large margin and furthermore shows improved or comparable
performance in the other tasks. Better convergence is a proven property of
optimal curricula [20]. In the lung cancer experiment, we showed that the
performance gain per training voxel is drastically increased compared to CPS.
This was observed for all different training configurations, differing in the
length of network training. The beneficial impact of our curriculum is based
on the ordering of patches by growing size, which we assume results in an
efficient, easy-to-hard order. Furthermore, PGPS also outperforms random
ordering (RPSS), validating the impact of the patch size ordering. We
hypothesize that with growing patch size, we ask the network subsequently more
difficult questions. In contrast, Li et al. [12] used only two different patch
sizes during training as a comparison method to their proposed curriculum.
Their results showed a performance drop of their patch-to-whole training over
CPS training. We assume that our approach benefits from maximal patch size
stages.

Moreover, our curriculum can utilize larger batch sizes within the same GPU
hardware restrictions. In our experiments, this resulted in better or
comparable performance for PGPS+ compared to classical CPS training, as higher
batch sizes generally can lead to better convergence and performance.

Besides the improved convergence, our proposed curriculum also drastically
reduced the training runtime and CO2 emissions to roughly half of the original
values for PGPS and PGPS+. This is a result of the cheaper operations due to
smaller input sizes. Minimal training runtime is defined by the total number
of seen voxels during training and depends highly on the technical
implementation.

We acknowledge that this is a proof of concept work, relying on the
hyperparameters of nnU-Net. More experiments are needed to explore the
different effects of using another minimal patch size, number of patch size
stages, batch size, normalization, and training scheduler. Furthermore, our
approach is not limited to fully convolutional networks but could also be
applied to transformers [3, 13] and hybrid model architectures [4, 5] as well
as object detection tasks.

Overall, the proposed Progressive Growing of Patch Size is a resource-
efficient curriculum strategy that drastically reduces the training runtime
and furthermore can also lead to improved convergence. The performance gain
per voxel is substantially higher than for classical constant patch size
training, supporting our hypothesis of better convergence. Improving the
efficiency of network training is particularly important considering the
growing carbon footprint of training deep learning models [17].

## References

  * [1] Antonelli, M., Reinke, A., Bakas, S., Farahani, K., Kopp-Schneider, A., Landman, B.A., Litjens, G., Menze, B., Ronneberger, O., Summers, R.M., et al.: The medical segmentation decathlon. Nature Communications 13(1), 4128 (2022). https://doi.org/10.1038/s41467-022-30695-9
  * [2] Bengio, Y., Louradour, J., Collobert, R., Weston, J.: Curriculum learning. In: Proceedings of the International Conference on Machine Learning. pp. 41–48. PMLR (2009)
  * [3] Dosovitskiy, A., Beyer, L., Kolesnikov, A., Weissenborn, D., Zhai, X., Unterthiner, T., Dehghani, M., Minderer, M., Heigold, G., Gelly, S., Uszkoreit, J., Houlsby, N.: An image is worth 16x16 words: Transformers for image recognition at scale. arXiv preprint arXiv:2010.11929 (2020)
  * [4] Hatamizadeh, A., Nath, V., Tang, Y., Yang, D., Roth, H.R., Xu, D.: Swin UNETR: Swin transformers for semantic segmentation of brain tumors in MRI images. In: Brainlesion Workshop at International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 272–284. Springer (2021)
  * [5] Hatamizadeh, A., Tang, Y., Nath, V., Yang, D., Myronenko, A., Landman, B., Roth, H.R., Xu, D.: UNETR: Transformers for 3D medical image segmentation. In: Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision. pp. 574–584 (2022)
  * [6] Havaei, M., Guizard, N., Chapados, N., Bengio, Y.: HeMIS: Hetero-modal image segmentation. In: Proceedings of International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 469–477. Springer (2016)
  * [7] Isensee, F., Jaeger, P.F., Kohl, S.A., Petersen, J., Maier-Hein, K.H.: nnU-Net: A self-configuring method for deep learning-based biomedical image segmentation. Nature Methods 18(2), 203–211 (2021)
  * [8] Isensee, F., Jäger, P.F., Kohl, S.A., Petersen, J., Maier-Hein, K.H.: Automated design of deep learning methods for biomedical image segmentation. arXiv preprint arXiv:1904.08128 (2019)
  * [9] Jesson, A., Guizard, N., Ghalehjegh, S.H., Goblot, D., Soudan, F., Chapados, N.: Cased: Curriculum adaptive sampling for extreme data imbalance. In: Proceedings of International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 639–646. Springer (2017)
  * [10] Jiménez-Sánchez, A., Mateus, D., Kirchhoff, S., Kirchhoff, C., Biberthaler, P., Navab, N., González Ballester, M.A., Piella, G.: Medical-based deep curriculum learning for improved fracture classification. In: Proceedings of International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 694–702. Springer (2019)
  * [11] Karras, T., Aila, T., Laine, S., Lehtinen, J.: Progressive growing of GANs for improved quality, stability, and variation. arXiv preprint arXiv:1710.10196 (2017)
  * [12] Li, H., Liu, X., Boumaraf, S., Liu, W., Gong, X., Ma, X.: A new three-stage curriculum learning approach for deep network based liver tumor segmentation. In: Proceedings of International Joint Conference on Neural Networks. pp. 1–6. IEEE (2020)
  * [13] Liu, Z., Lin, Y., Cao, Y., Hu, H., Wei, Y., Zhang, Z., Lin, S., Guo, B.: Swin transformer: Hierarchical vision transformer using shifted windows. In: Proceedings of the IEEE/CVF International Conference on Computer Vision. pp. 10012–10022 (2021)
  * [14] Platanios, E.A., Stretcu, O., Neubig, G., Poczos, B., Mitchell, T.M.: Competence-based curriculum learning for neural machine translation. arXiv preprint arXiv:1903.09848 (2019)
  * [15] Ronneberger, O., Fischer, P., Brox, T.: U-Net: Convolutional networks for biomedical image segmentation. In: Proceedings of International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 234–241. Springer (2015)
  * [16] Schmidt, V., Goyal, K., Joshi, A., Feld, B., Conell, L., Laskaris, N., Blank, D., Wilson, J., Friedler, S., Luccioni, S.: CodeCarbon: Estimate and track carbon emissions from machine learning computing (2021). https://doi.org/10.5281/zenodo.4658424, v2.3.4
  * [17] Selvan, R., Bhagwat, N., Wolff Anthony, L.F., Kanding, B., Dam, E.B.: Carbon footprint of selecting and training deep learning models for medical image analysis. In: Proceedings of International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 506–516. Springer (2022)
  * [18] Spitkovsky, V.I., Alshawi, H., Jurafsky, D.: From baby steps to leapfrog: How “less is more” in unsupervised dependency parsing. In: Proceedings of Annual Conference of the North American Chapter of the Association for Computational Linguistics. pp. 751–759. Association for Computational Linguistics (2010)
  * [19] Wei, J., Suriawinata, A., Ren, B., Liu, X., Lisovsky, M., Vaickus, L., Brown, C., Baker, M., Nasir-Moin, M., Tomita, N., Torresani, L., Wei, J., Hassanpour, S.: Learn like a pathologist: Curriculum learning by annotator agreement for histopathology image classification. In: Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision. pp. 2473–2483 (2021)
  * [20] Weinshall, D., Cohen, G., Amir, D.: Curriculum learning by transfer learning: Theory and experiments with deep networks. In: Proceedings of International Conference on Machine Learning. pp. 5238–5246. PMLR (2018)
  * [21] Zaremba, W., Sutskever, I.: Learning to execute. arXiv preprint arXiv:1410.4615 (2014)
  * [22] Zhao, J., Dai, L., Zhang, M., Yu, F., Li, M., Li, H., Wang, W., Zhang, L.: PGU-net+: Progressive growing of U-net+ for automated cervical nuclei segmentation. In: Multiscale Multimodal Medical Imaging Workshop at International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 51–58. Springer (2020)

Supplemental Table 2: Performance of Constant Patch Size Training Scheme (CPS) and Progressive Growing of Patch Size Training Scheme (PGPS) and Random Patch Size Sampling Training Scheme (RPSS) on Task06 Lung Cancer of Medical Segmentation Decathlon for different numbers of training iterations. (Iterations Per Epoch: Number of training iterations per epoch normalized to default training iterations of nnU-Net; Dice Score: Evaluated in 5-fold Cross-Validation; Runtime: Training runtime per fold on one NVIDIA A100-SXM4-80GB GPU; Voxels Shown: Total number of voxels iterated during training normalized to CPS training) Iterations Per Epoch | Training Scheme | Dice Score | Runtime | Voxels Shown
---|---|---|---|---
10% | CPS | 0.6073 | 2.49 h | 10.0 %
| RPSS | 0.5971 | 1.72 h |  $\sim$ 3.5 %
| PGPS | 0.6608 | 1.69 h | 3.5 %
25% | CPS | 0.6435 | 4.22 h | 25.0 %
| RPSS | 0.6386 | 2.37 h |  $\sim$ 8.8 %
| PGPS | 0.6753 | 2.35 h | 8.8 %
50% | CPS | 0.6833 | 7.40 h | 50.0 %
| RPSS | 0.7005 | 3.52 h |  $\sim$ 17.7 %
| PGPS | 0.7009 | 3.43 h | 17.7 %
100% | CPS | 0.7211 | 13.55 h | 100.0 %
| RPSS | 0.7173 | 5.69 h |  $\sim$ 35.4 %
| PGPS | 0.7263 | 5.70 h | 35.4 %
Supplemental Table 3: All training input tensor sizes for nnU-Net model
training with Progressive Growing of Patch Size Curriculum on the Medical
Segmentation Decathlon tasks. The input is given as batch size * width *
height * depth. These are the values used for PGPS+. The raw PGPS uses a
constant batch size for all patch size stages, which is here given for the
final patch size. The maximum patch size during model training is also used
for inference.

|  Brain |  Heart |  Liver |  Hippocampus |  Prostate
---|---|---|---|---|---
1\. Stage |  24*64*32*32 |  24*32*32*32 |  24*64*32*32 |  24*16*8*8 |  24*8*64*64
2\. Stage |  12*64*64*32 |  12*32*64*32 |  12*64*64*32 |  12*16*16*8 |  12*8*128*64
3\. Stage |  6*64*64*64 |  6*32*64*64 |  6*64*64*64 |  9*16*16*16 |  6*8*128*128
4\. Stage |  4*96*64*64 |  4*48*64*64 |  4*96*64*64 |  9*24*16*16 |  4*12*128*128
5\. Stage |  3*96*96*64 |  3*48*96*64 |  3*96*96*64 |  9*24*24*16 |  3*12*192*128
6\. Stage |  2*96*96*96 |  2*48*96*96 |  2*96*96*96 |  9*24*24*24 |  2*12*192*192
7\. Stage |  2*128*96*96 |  2*64*96*96 |  2*128*96*96 |  9*32*24*24 |  2*16*192*192
8\. Stage |  2*128*128*96 |  2*64*128*96 |  2*128*128*96 |  9*32*32*24 |  2*16*256*192
9\. Stage |  2*128*128*128 |  2*64*128*128 |  2*128*128*128 |  9*32*32*32 |  2*16*256*256
10\. Stage |  - |  2*80*128*128 |  - |  9*40*32*32 |  2*20*256*256
11\. Stage |  - |  2*80*160*128 |  - |  9*40*40*32 |  2*20*320*256
12\. Stage |  - |  2*80*160*160 |  - |  9*40*40*40 |  -
13\. Stage |  - |  2*80*192*160 |  - |  9*40*48*40 |  -
14\. Stage |  - |  - |  - |  9*40*56*40 |  -
15\. Stage |  - |  - |  - |  - |  -

|  Lung |  Pancreas |  Hepatic Vessel |  Spleen |  Colon
---|---|---|---|---|---
1\. Stage |  24*32*32*32 |  24*16*32*32 |  24*32*32*32 |  24*32*32*32 |  24*16*32*32
2\. Stage |  12*32*64*32 |  12*16*64*32 |  12*32*64*32 |  12*32*64*32 |  12*16*64*32
3\. Stage |  6*32*64*64 |  6*16*64*64 |  6*32*64*64 |  6*32*64*64 |  6*16*64*64
4\. Stage |  4*48*64*64 |  4*24*64*64 |  4*48*64*64 |  4*48*64*64 |  4*24*64*64
5\. Stage |  3*48*96*64 |  3*24*96*64 |  3*48*96*64 |  3*48*96*64 |  3*24*96*64
6\. Stage |  2*48*96*96 |  2*24*96*96 |  2*48*96*96 |  2*48*96*96 |  2*24*96*96
7\. Stage |  2*64*96*96 |  2*32*96*96 |  2*64*96*96 |  2*64*96*96 |  2*32*96*96
8\. Stage |  2*64*128*96 |  2*32*128*96 |  2*64*128*96 |  2*64*128*96 |  2*32*128*96
9\. Stage |  2*64*128*128 |  2*32*128*128 |  2*64*128*128 |  2*64*128*128 |  2*32*128*128
10\. Stage |  2*80*128*128 |  2*40*128*128 |  2*64*160*128 |  2*64*160*128 |  2*40*128*128
11\. Stage |  2*80*160*128 |  2*40*160*128 |  2*64*160*160 |  2*64*160*160 |  2*40*160*128
12\. Stage |  2*80*160*160 |  2*40*160*160 |  2*64*192*160 |  2*64*192*160 |  2*40*160*160
13\. Stage |  2*80*192*160 |  2*40*192*160 |  2*64*192*192 |  - |  2*48*160*160
14\. Stage |  - |  2*40*192*192 |  - |  - |  2*48*192*160
15\. Stage |  - |  2*40*224*192 |  - |  - |  2*56*192*160


* Conti Nibali et al. 2023 Conti Nibali, V.; Maiti, S.; Saija, F.; Heyden, M.; Cassone, G. Electric-field induced entropic effects in liquid water. _The Journal of Chemical Physics_ 2023, _158_ , 184501
  * Ojha and Kühne 2023 Ojha, D.; Kühne, T. D. Vibrational dynamics of liquid water in an external electric field. _Phys. Chem. Chem. Phys._ 2023, _25_ , 13442–13451
  * Cassone and Martelli 2024 Cassone, G.; Martelli, F. Electrofreezing of liquid water at ambient conditions. _Nature Communications_ 2024, _15_ , 1856
  * Persson and Halle 2015 Persson, F.; Halle, B. How amide hydrogens exchange in native proteins. _Proceedings of the National Academy of Sciences_ 2015, _112_ , 10383–10388
  * Zhang and Jiang 2023 Zhang, Y.; Jiang, B. Universal machine learning for the response of atomistic systems to external fields. _Nature Communications_ 2023, _14_ , 6424
  * Dral and Barbatti 2021 Dral, P. O.; Barbatti, M. Molecular excited states through a machine learning lens. _Nature Reviews Chemistry_ 2021, _5_ , 388–405
  * Molinero and Moore 2009 Molinero, V.; Moore, E. B. Water Modeled As an Intermediate Element between Carbon and Silicon. _The Journal of Physical Chemistry B_ 2009, _113_ , 4008–4016
  * Moore and Molinero 2011 Moore, E. B.; Molinero, V. Structural transformation in supercooled water controls the crystallization rate of ice. _Nature_ 2011, _479_ , 506–508
  * Holten et al. 2013 Holten, V.; Limmer, D. T.; Molinero, V.; Anisimov, M. A. Nature of the anomalies in the supercooled liquid state of the mW model of water. _The Journal of Chemical Physics_ 2013, _138_ , 174501
  * Xu and Molinero 2011 Xu, L.; Molinero, V. Is There a Liquid–Liquid Transition in Confined Water? _The Journal of Physical Chemistry B_ 2011, _115_ , 14210–14216
  * Moore et al. 2010 Moore, E. B.; de la Llave, E.; Welke, K.; Scherlis, D. A.; Molinero, V. Freezing, melting and structure of ice in a hydrophilic nanopore. _Phys. Chem. Chem. Phys._ 2010, _12_ , 4124–4134
  * Factorovich et al. 2015 Factorovich, M. H.; Molinero, V.; Scherlis, D. A. Hydrogen-Bond Heterogeneity Boosts Hydrophobicity of Solid Interfaces. _Journal of the American Chemical Society_ 2015, _137_ , 10618–10623
  * Factorovich et al. 2014 Factorovich, M. H.; Molinero, V.; Scherlis, D. A. A simple grand canonical approach to compute the vapor pressure of bulk and finite size systems. _The Journal of Chemical Physics_ 2014, _140_ , 064111
  * Gadea et al. 2020 Gadea, E. D.; Perez Sirkin, Y. A.; Molinero, V.; Scherlis, D. A. Electrochemically Generated Nanobubbles: Invariance of the Current with Respect to Electrode Size and Potential. _The Journal of Physical Chemistry Letters_ 2020, _11_ , 6573–6579
  * Jacobson et al. 2010 Jacobson, L. C.; Hujo, W.; Molinero, V. Amorphous Precursors in the Nucleation of Clathrate Hydrates. _Journal of the American Chemical Society_ 2010, _132_ , 11806–11811
  * DeMille and Molinero 2009 DeMille, R. C.; Molinero, V. Coarse-grained ions without charges: Reproducing the solvation structure of NaCl in water using short-ranged potentials. _The Journal of Chemical Physics_ 2009, _131_ , 034107
  * Perez Sirkin et al. 2016 Perez Sirkin, Y. A.; Factorovich, M. H.; Molinero, V.; Scherlis, D. A. Vapor Pressure of Aqueous Solutions of Electrolytes Reproduced with Coarse-Grained Models without Electrostatics. _Journal of Chemical Theory and Computation_ 2016, _12_ , 2942–2949
  * Lu et al. 2019 Lu, J.; Barnett, A.; Molinero, V. Effect of Polymer Architecture on the Nanophase Segregation, Ionic Conductivity, and Electro-Osmotic Drag of Anion Exchange Membranes. _The Journal of Physical Chemistry C_ 2019, _123_ , 8717–8726
  * Mochizuki and Molinero 2018 Mochizuki, K.; Molinero, V. Antifreeze Glycoproteins Bind Reversibly to Ice via Hydrophobic Groups. _Journal of the American Chemical Society_ 2018, _140_ , 4803–4811
  * Chan et al. 2019 Chan, H.; Cherukara, M. J.; Narayanan, B.; Loeffler, T. D.; Benmore, C.; Gray, S. K.; Sankaranarayanan, S. K. R. S. Machine learning coarse grained models for water. _Nat. Commun._ 2019, _10_


#  tsdownsample: high-performance time series downsampling for scalable
visualization

Jeroen Van Der Donckt  Jonas Van Der Donckt  Sofie Van Hoecke  IDLab, Ghent
University - imec, Technologiepark Zwijnaarde 126, 9052 Zwijnaarde, Belgium

###### Abstract

Interactive line chart visualizations greatly enhance the effective
exploration of large time series. Although downsampling has emerged as a well-
established approach to enable efficient interactive visualization of large
datasets, it is not an inherent feature in most visualization tools.
Furthermore, there is no library offering a convenient interface for high-
performance implementations of prominent downsampling algorithms. To address
these shortcomings, we present tsdownsample, an open-source Python package
specifically designed for CPU-based, in-memory time series downsampling. Our
library focuses on performance and convenient integration, offering optimized
implementations of leading downsampling algorithms. We achieve this
optimization by leveraging low-level SIMD instructions and multithreading
capabilities in Rust. In particular, SIMD instructions were employed to
optimize the argmin and argmax operations. This SIMD optimization, along with
some algorithmic tricks, proved crucial in enhancing the performance of
various downsampling algorithms. We evaluate the performance of tsdownsample
and demonstrate its interoperability with an established visualization
framework. Our performance benchmarks indicate that the algorithmic runtime of
tsdownsample approximates the CPU’s memory bandwidth. This work marks a
significant advancement in bringing high-performance time series downsampling
to the Python ecosystem, enabling scalable visualization. The open-source code
can be found at https://github.com/predict-idlab/tsdownsample

###### keywords:

time series , visualization , downsampling , Rust , Python , SIMD

††journal: SoftwareX

##  1 Introduction

Time series are ubiquitous in many domains, such as healthcare, finance, and
manufacturing. This complex data modality can be challenging to comprehend
through summary statistics alone, making visualizations a crucial tool for
gaining insights, with line charts proving particularly effective for most
tasks [1]. By following the “overview first, zoom and filter, then details on
demand” paradigm [2], interactive line chart visualizations allow users to
quickly and easily understand the data and identify patterns and trends [3].

Many real-world time series datasets are extremely large, encompassing
millions or even billions of data points. As a result, there is a pressing
need for scalable visualization techniques that are capable of effectively
handling such datasets [4, 5]. One approach to realize scalable visualization
is through utilizing data aggregation techniques such as downsampling, which
reduces the number of data points in a time series while preserving its
overall shape [6, 7, 1]. Downsampling enables faster rendering and more
responsive interactions, allowing users to explore large datasets more
effectively [8, 9]. Downsampling algorithms find wide adoption in the time
series database domain, with Uber integrating a downsampling function in their
M3 metrics platform [10], and TimeScaleDB offering downsampling as a server-
side hyperfunction [11].

However, when dealing with very large datasets (billions of data points), the
downsampling process itself can become a bottleneck [12]. This is especially
the case in the context of interactive visualizations, which require fast
downsampling to minimize latency when interacting with the graph, such as
zooming and panning [9]. Moreover, the authors observe that, at the time of
writing, there is no Python library offering high-performance implementations
of multiple time series downsampling algorithms.

Recognizing this challenge, we introduce tsdownsample, an open-source Python
toolkit designed for in-memory, CPU-based, time series downsampling, focusing
on performance and integrability. tsdownsample provides optimized CPU
implementations of the most prominent downsampling algorithms, i.e., EveryNth,
MinMax, M4 [7], LTTB (Largest-Triangle-Three-Buckts) [6], and MinMaxLTTB [12].
The algorithms are implemented in Rust, a system programming language known
for its performance and memory safety. The Rust code leverages SIMD (Single
Instruction Multiple Data) instructions together with some algorithmic tricks
and (optionally) multithreading to achieve exceptional performance and
scalability. A core component of tsdownsample is our ArgMinMax Rust library,
which provides SIMD accelerated argmin and argmax functionality. Optimizing
these operations for various CPU architectures proved to be crucial, as they
form the inner loop of most downsampling algorithms [13]. tsdownsample is
distributed as a Python toolkit by publishing the cross-compiled Python
bindings for the underlying Rust code for a wide range of operating systems
and CPU architectures.

In summary, this paper contributes tsdownsample, a high-performance library
optimized for CPU that provides downsampling for scalable time series
visualization. This library’s integrability is demonstrated through its
adoption as the downsampling solution in a time series visualization library,
which has over 1.5 million installations at the time of writing.

The remainder of this paper is structured as follows. In Section 2, we present
a description of the software. We then provide an illustrative example in
Section 3, where we also highlight its integration in an existing toolkit.
Finally, we evaluate the performance of tsdownsample in Section 4.

##  2 Software description

tsdownsample is a Python package that utilizes Rust to provide CPU-optimized
implementations of downsampling algorithms for time series visualization. To
facilitate seamless installation and usage, Python bindings for the underlying
Rust code are cross-compiled for various operating systems and CPU
architectures, which are distributed as a PyPi package. Installing
tsdownsample is simple and can be done via pip by running the following
command: pip install tsdownsample.

The following subsections detail the performance and user-interface aspect of
tsdownsample. We will first describe how optimizing argmin and argmax
operations in Rust proved crucial for achieving high-performance downsampling.
Then, we will discuss the convenient interface that this library provides for
all implemented downsampling algorithms.

###  2.1 ArgMinMax

Given the significance of vertical extrema for ensuring the visual
representativeness of time series downsampling [13], it was imperative to
optimize the argmin and argmax operations. These operations play a vital role
in the inner loop of the MinMax, M4 [7], and MinMaxLTTB [12] algorithms. As
such, we developed the ArgMinMax Rust library (also referred to as a crate),
which provides a highly efficient and overflow-free implementation of the
argmin and argmax operations. These operations return the indices of the
minimum and maximum values of an array. Note that the argmin and argmax values
are extracted simultaneously within a single pass over the data, as this is
mainly a memory-bound task.

The ArgMinMax crate includes SIMD-optimized implementations of argmin and
argmax for SSE, AVX(2), AVX512, and NEON, and includes runtime CPU feature
detection to select the optimal (supported) SIMD implementation for the
current CPU. SIMD instructions allow the CPU to perform the same operation on
multiple data points simultaneously, providing a significant boost in
performance for certain types of operations. The library is SIMD-optimized for
a wide range of CPU architectures; x86, x86_64, arm(v7), and aarch64. In
addition, the ArgMinMax crate supports a wide range of data types (f16, f32,
f64, i8, i16, i32, i64, u8, u16, u32, and u64). We further guarantee the
library to be memory-efficient, as it operates on a memory view (i.e., a
slice) of the data rather than copying it. The SIMD algorithm is also
branchless, ensuring that the runtime is independent of the quality of the
branch predictor, making the best-case runtime the same as the worst-case
runtime.

####  2.1.1 argminmax SIMD algorithm

In code snippet 1, we present the inner loop of the SIMD argmin and argmax
algorithm. This algorithm extracts both the argmin and argmax value in a
single pass over the data. To do so, we utilize four accumulating SIMD vectors
(also referred to as registers). Two of these registers maintain the lowest
and highest values encountered while iterating over the data in chunks of size
LANE_SIZE. At the end of the iteration, these vectors contain the maximum and
minimum values at each position within all seen LANE_SIZE chunks. As a final
step, after iterating over all the chunks, the algorithm extracts the minimum
and maximum values along with their respective indices from the SIMD vectors
(i.e., the horizontal operations).

⬇

arr_ptr = arr.as_ptr(); // Array pointer we will increment in the loop

new_index = INITIAL_INDEX; // Index we will increment in the loop

// Initialization of the accumulating SIMD vectors

// Note that for adequate float (NaN) handling this initialization should be
different

index_low = INITIAL_INDEX;

values_low = _mm_loadu(arr_ptr);

index_high = INITIAL_INDEX;

values_high = _mm_loadu(arr_ptr);

for _ in 0..arr.len() / LANE_SIZE - 1 { // Iterate over the array in LANE_SIZE
chunks

// Increment the index

new_index = _mm_add(new_index, INDEX_INCREMENT);

// Load the next chunk of data

arr_ptr = arr_ptr.add(LANE_SIZE);

new_values = _mm_loadu(arr_ptr);

// Update the lowest values and index

mask_low = _mm_cmplt(new_values, values_low);

values_low = _mm_blendv(values_low, new_values, mask_low);

index_low = _mm_blendv(index_low, new_index, mask_low);

// Update the highest values and index

mask_high = _mm_cmpgt(new_values, values_high);

values_high = _mm_blendv(values_high, new_values, mask_high);

index_high = _mm_blendv(index_high, new_index, mask_high);

}

// Get the min/max index and corresponding value from the SIMD vectors

(min_index, min_value) = _horiz_min(index_low, values_low);

(max_index, max_value) = _horiz_max(index_high, values_high);

Figure 1: The core (inner loop) of the SIMD argminmax algorithm.

The pseudocode in snippet 1 closely resembles the Rust code of the ArgMinMax
package. It is worth noting that the SIMD instructions, such as _mm_loadu,
_mm_cmplt, _mm_cmpgt, and _mm_blendv, are generic function names that need to
be associated with the corresponding CPU instructions of the various
architectures. To achieve this in Rust, we utilized a trait that defines these
generic SIMD instructions as functions, similar to C++ templates. For each
supported CPU architecture and data type combination in the ArgMinMax package,
we have implemented a concrete version of this trait111Since the lane-size is
also influenced by this combination, this trait includes the lane size
parameter as well..

When the length of the array exceeds the maximum value that the index vector’s
underlying data type can represent, an (index) overflow will occur. It is this
overflow challenge that makes the argmin and argmax operations a much harder
problem to SIMD-optimize compared to the min and max operations. As a result,
compiling a scalar implementation to vectorized instructions is not trivial or
even impossible. As such, it was necessary to manually write the algorithm
using SIMD instructions, rather than relying on the compiler for optimized
compilation. Notably, the polars library, an exceptionally fast DataFrame
library and in-memory query engine that currently exceeds 0.7M monthly
installations, adopted our ArgMinMax crate to provide more optimized argmin
and argmax operations [14]. For more details on how to implement an overflow-
free solution (which requires an additional outer loop), please refer to the
open-source code repository available at github.com/predict-idlab/argminmax.

####  2.1.2 Optimized implementation for f16 and uints

In contrast to other data types, most modern CPUs (x86) do not have hardware
support for the float16 (f16) data type222With float16 we refer to IEEE
754-2008 standard binary16, also known as half floating point type [15].. As a
result, programming languages typically support f16 either by upcasting to f32
or by using a software implementation. Both approaches come at the cost of
considerable overhead.

Instead of applying one of these two approaches, we convert f16 to an ordinal
mapping of i16 (which we refer to as i16ord). This allows us to efficiently
support f16 data types in the ArgMinMax crate and in the tsdownsample package
as a whole. The mapping preserves the ordinality of the f16 data, as
illustrated in Figure 3, allowing the use of fast built-in i16 (SIMD)
instructions for comparison333Note that this mapping only makes sense when you
are solely interested in comparing values.. Moreover, the transformation is
symmetric, meaning that we can transform the outcome back to f16 (by using the
same mapping function) without needing a lookup table, as illustrated in
Figure 3. Furthermore, as the transformation only performs binary (bitwise)
operations, the overhead is limited, thus making it convenient to implement
using SIMD instructions.

The ordinal transformation is performed as follows444To apply this
transformation to a f16 value, we first transmute the f16 value to i16.:

ord_transform(v: i16) = ((v >> 15) & 0x7FFF) $\oplus$ v)

Figure 2: f16 $\rightarrow$ i16

Figure 3: f16 $\rightarrow$ i16 $\rightarrow$ f16

A final implementation remark is that uint data types (u8, u16, u32, u64) lack
SSE, AVX(2), nor AVX512 SIMD instructions. This data type is supported by
performing a similar, albeit much simpler, ordinal mapping to integer data
type. This mapping leverages the two’s complement of signed integers. In
particular, the transformation first transmutes the uint as int and then XORs
the value with the smallest (i.e., the largest negative) int value.

###  2.2 tsdownsample

tsdownsample further builds upon the optimizations of the ArgMinMax crate. In
particular, the MinMax, M4, and MinMaxLTTB algorithms directly rely on the
argmingmax algorithm in their inner loop (i.e., for each bin). Since these
algorithms operate on local heuristics within each bin, they can be easily
parallelized in Rust to leverage the processing power of modern multicore CPUs
[16]. To achieve this, we implemented a multithreaded bin index generator
using a search sorted approach, thereby enhancing cache hits through chunked
execution. Remark that multithreading is not possible in Python due to the
Global Interpreter Lock (GIL) which prevents multiple threads from executing
Python bytecodes simultaneously [17].

Unit testing is conducted to ensure the correctness of the supported
downsampling algorithms. Specifically, we verify the consistency of the
downsampled (i.e., selected) data points across various data types,
downsamplers, and multithreading configurations. Additionally, we compare the
Rust implementations to a reference Python implementation, which (although
being considerably slower) serves as a benchmark for correctness. Furthermore,
we ensure that passing an equally sampled x yields the same output as not
specifying the index. Noteworthy, to capture potential performance regressions
when updating the code base, we added performance monitoring to the CI/CD
workflow.

###  2.3 Downsampling interface

tsdownsample aims to provide a convenient interface for the supported
downsampling algorithms. Users can interact with these algorithms through
Python classes, which act as a thin wrapper around the underlying Rust
bindings. These classes abstract the dispatching of data type specific
function calls, making the interface easy to use. All classes implement the
downsample method, which has the following signature:

downsample([x], y, n_out, **kwargs) -> ndarray[uint64]

This signature first accepts two positional arguments that represent the input
values555This signature design aligns with the convention used in the
matplotlib.pyplot.plot method [18].. The first positional argument x is
optional and represents the index of the time series values (y). If not
provided, it is assumed that the time series values are equally sampled
without any gaps. The second positional argument y is mandatory and
corresponds to the input time series values. The n_out argument is a mandatory
keyword argument that defines the number of output values666It is important to
note that if there are gaps in the time series (index), fewer than n_out
indices may be returned, as no data points can be selected for empty bins. In
addition to these arguments, optional keyword arguments can be passed via
**kwargs. These additional arguments provide increased flexibility, including
options such as the parallel argument, a boolean that enables multi-threading
when set to True. By default, the parallel option is set to False. Lastly, the
downsample method returns a numpy array containing unsigned 64-bit integers,
representing the indices of the downsampled (i.e., selected) values.

##  3 Illustrative example

Listing 4 provides an illustration of how tsdownsample can be utilized. The
code snippet begins by importing the MinMaxLTTBDownsampler class from the
tsdownsample package, along with the numpy library. numpy is utilized to
generate a random time series dataset consisting of 10 million points. Next,
an instance of the MinMaxLTTBDownsampler class is constructed, which is used
to downsample the aforementioned time series data to 1000 points. This is
accomplished by utilizing the downsample method, whose interface is detailed
in Section 2.3. The resulting indices of the downsampled time series data are
stored in the s_ds variable. Subsequently, these selected indices can be used
to retrieve a representative subset of the original time series data,
facilitating efficient visualization.

⬇

from tsdownsample import MinMaxLTTBDownsampler

import numpy as np

# Create a time series

y = np.random.randn(10_000_000)

# Downsample to 1000 points

# ———————————-

# Get the selected indices when downsampling y to 1000 points

s_ds = MinMaxLTTBDownsampler().downsample(y, n_out=1000)

Figure 4: Downsampling a random array with MinMaxLTTB.

###  3.1 Integration in plotly-resampler

tsdownsample has been integrated as downsampling back end in the plotly-
resampler visualization tool since version v0.9 [8]. At the time of writing,
plotly-resampler has over 1.5 million installations. One longstanding
challenge with the plotly-resampler library was the need for users to compile
downsampling C code locally during the installation process. This requirement
often led to complications, as it necessitates users to have the correct
Python headers, ensure compatibility of the numpy version with the utilized C
API, and have an appropriate C compiler installed. However, by adopting
tsdownsample, these issues have been effectively addressed. tsdownsample has a
precompiled binary available for multiple platforms, eliminating the need for
users to compile the underlying code themselves, while ensuring optimal
performance. This integration has not only resolved these compilation-related
obstacles but has also resulted in remarkable speed improvements, ranging from
3 to 30 times faster performance.

##  4 Performance

We analyzed the performance of tsdownsample for a range of data types and
algorithms, as presented in Table 1. To do so, we created an array of random
values for the data types under consideration and measured the time required
to downsample the respective array to 2,000 values (i.e., n_out=2,000) using
Python’s timeit module. A reproducible notebook containing the benchmark code
can be found here: https://github.com/predict-
idlab/tsdownsample/blob/main/notebooks/benches.ipynb.

The benchmarks were executed on a server with an Intel Xeon E5-2650 v2 (32) @
3.40GHz CPU and SAMSUNG M393B1G73QH0-CMA DDR3 1600MT/s RAM, running on the
Ubuntu 18.04.6 LTS x86_64 operating system. Other running processes were
limited to a minimum.

Table 1 displays the median time measurements for the five downsampling
algorithms available in tsdownsample. These measurements are provided for all
supported data types and varying numbers of data points. Among the algorithms,
EveryNth exhibits a constant execution time of approximately 0.02 ms,
regardless of the length of the input data. For the M4, MinMax, and MinMaxLTTB
algorithms, we observe two main trends. Firstly, there is a linear or
sublinear increase (approximately 10x or less) in runtime when transitioning
from 10 million to 100 million and then to 1 billion data points. Secondly,
these three algorithms have similar runtimes for the same data type, which is
most noticeable in the 1 billion data point rows. Both these trends hold true
for the sequential and parallel executions. In the case of the LTTB algorithm,
we again notice a linear scaling pattern as the data length increases. It is
important to note that the runtime of LTTB is significantly slower, up to two
orders of magnitude, compared to MinMaxLTTB, especially when dealing with
larger datasets (e.g., 1 billion uint8 data points). This slower runtime can
be attributed to LTTB requiring much more computationally expensive
calculations [6], which is largely mitigated in MinMaxLTTB [12].

| Algorithm | EveryNth | M4 | MinMax | MinMaxLTTB | LTTB
---|---|---|---|---|---|---
| Parallel | False | False | True | False | True | False | True | False
dtype | N |  |  |  |  |  |  |  |
float16 | 1,000,000 | 0.02 | 0.47 | 0.43 | 0.43 | 0.55 | 1.12 | 0.80 | 6.60
10,000,000 | 0.01 | 2.94 | 0.59 | 2.31 | 0.56 | 4.13 | 0.89 | 59.95
100,000,000 | 0.03 | 24.93 | 4.94 | 34.49 | 4.97 | 25.89 | 5.25 | 575.40
1,000,000,000 | 0.01 | 255.13 | 44.61 | 250.30 | 44.92 | 262.83 | 45.06 | 5614.58
float32 | 1,000,000 | 0.01 | 0.41 | 0.27 | 0.46 | 0.23 | 1.17 | 0.56 | 2.41
10,000,000 | 0.03 | 4.02 | 0.89 | 3.56 | 0.94 | 6.29 | 1.15 | 18.39
100,000,000 | 0.02 | 41.27 | 9.25 | 33.40 | 9.28 | 40.40 | 9.56 | 173.72
1,000,000,000 | 0.02 | 407.50 | 88.64 | 338.46 | 88.75 | 398.32 | 89.17 | 1750.17
float64 | 1,000,000 | 0.02 | 0.73 | 0.33 | 0.75 | 0.34 | 1.52 | 0.53 | 2.33
10,000,000 | 0.02 | 8.60 | 2.04 | 8.80 | 2.18 | 10.51 | 2.30 | 18.25
100,000,000 | 0.02 | 85.43 | 18.06 | 83.51 | 18.09 | 86.07 | 18.37 | 196.97
1,000,000,000 | 0.01 | 832.17 | 176.74 | 663.86 | 176.82 | 828.15 | 177.26 | 1804.99
int8 | 1,000,000 | 0.03 | 0.38 | 0.43 | 0.46 | 0.45 | 0.93 | 0.74 | 3.40
10,000,000 | 0.02 | 1.98 | 0.51 | 1.50 | 0.52 | 2.64 | 0.77 | 24.39
100,000,000 | 0.01 | 18.08 | 2.69 | 14.49 | 2.77 | 17.58 | 3.03 | 237.54
1,000,000,000 | 0.01 | 151.21 | 22.63 | 142.52 | 22.61 | 155.97 | 22.83 | 2374.65
int16 | 1,000,000 | 0.02 | 0.38 | 0.35 | 0.41 | 0.37 | 0.86 | 0.64 | 3.34
10,000,000 | 0.01 | 3.78 | 0.46 | 2.01 | 0.50 | 3.41 | 0.74 | 28.64
100,000,000 | 0.02 | 30.17 | 4.95 | 23.65 | 5.07 | 23.87 | 5.20 | 253.98
1,000,000,000 | 0.02 | 230.80 | 44.80 | 229.14 | 44.80 | 232.49 | 45.00 | 2417.14
int32 | 1,000,000 | 0.03 | 0.59 | 0.38 | 0.64 | 0.40 | 1.28 | 0.63 | 3.19
10,000,000 | 0.02 | 5.36 | 0.92 | 4.25 | 1.01 | 5.05 | 1.21 | 23.79
100,000,000 | 0.01 | 44.57 | 9.61 | 45.47 | 9.31 | 58.25 | 9.78 | 227.60
1,000,000,000 | 0.01 | 452.85 | 88.67 | 390.95 | 88.72 | 470.03 | 89.00 | 2297.82
int64 | 1,000,000 | 0.03 | 1.52 | 0.45 | 1.55 | 0.57 | 1.86 | 0.70 | 3.25
10,000,000 | 0.01 | 13.33 | 2.20 | 11.15 | 2.28 | 11.55 | 2.50 | 24.63
100,000,000 | 0.03 | 116.66 | 18.50 | 139.79 | 18.17 | 134.90 | 19.30 | 270.74
1,000,000,000 | 0.03 | 1154.85 | 177.20 | 1114.58 | 177.47 | 1171.36 | 177.52 | 2408.60
uint8 | 1,000,000 | 0.02 | 0.39 | 0.44 | 0.47 | 0.48 | 0.91 | 0.74 | 3.42
10,000,000 | 0.02 | 2.08 | 0.52 | 1.61 | 0.53 | 2.65 | 0.77 | 24.52
100,000,000 | 0.02 | 15.97 | 2.70 | 15.01 | 2.78 | 18.86 | 2.96 | 243.15
1,000,000,000 | 0.01 | 155.68 | 22.73 | 147.67 | 22.62 | 160.77 | 22.79 | 2384.09
uint16 | 1,000,000 | 0.02 | 0.52 | 0.34 | 0.37 | 0.35 | 0.81 | 0.55 | 3.29
10,000,000 | 0.01 | 2.56 | 0.45 | 2.02 | 0.45 | 3.22 | 0.62 | 24.86
100,000,000 | 0.01 | 24.31 | 4.88 | 20.95 | 4.82 | 24.79 | 4.89 | 317.30
1,000,000,000 | 0.01 | 239.21 | 44.52 | 202.76 | 44.70 | 238.73 | 44.74 | 2392.48
uint32 | 1,000,000 | 0.02 | 0.58 | 0.34 | 0.62 | 0.36 | 1.15 | 0.57 | 3.30
10,000,000 | 0.02 | 5.98 | 0.92 | 4.33 | 0.93 | 8.88 | 1.16 | 25.02
100,000,000 | 0.01 | 47.35 | 9.31 | 40.29 | 9.27 | 48.17 | 9.56 | 240.22
1,000,000,000 | 0.01 | 468.04 | 88.72 | 399.56 | 88.91 | 475.70 | 88.91 | 2402.57
uint64 | 1,000,000 | 0.01 | 1.85 | 0.40 | 1.87 | 0.42 | 2.14 | 0.61 | 3.77
10,000,000 | 0.01 | 14.69 | 2.11 | 13.90 | 2.15 | 14.17 | 2.42 | 35.38
100,000,000 | 0.02 | 140.70 | 18.19 | 137.31 | 18.21 | 140.87 | 18.36 | 276.40
1,000,000,000 | 0.01 | 1403.84 | 176.99 | 1380.90 | 177.19 | 1402.95 | 177.25 | 2817.33

Table 1: Downsampling time (in ms) for the algorithms offered by tsdownsample.
The experiment parameters are described in the first two columns: dtype
indicates the data type, and N denotes the number of data points. Note that
LLTB cannot be parallelized, since this algorithm requires a sequential
iteration over the bins [12]. Figure 5: Downsampling time comparison for
integer data types for the MinMax, M4, and MinMaxLTTB algorithms (shown in the
three subplots) provided by tsdownsample. The y-axis represents the
downsampling time in milliseconds on a logarithmic scale, while the x-axis
indicates the number of data points. The figure includes (nearly
imperceptible) whiskers to denote the standard deviation of the measurements
for the collected data points.

Figure 5 provides additional insights to complement the findings of Table 1.
It illustrates the relationship between downsampling time (y-axis) and the
number of data points (x-axis) for the three algorithms that utilize the
optimizations of ArgMinMax. Given that the y-axis is logarithmic scale, the
logarithmic trend that we observe for all integer data types in every subplot,
confirms the earlier mentioned observation that the implementation scales
linearly with the number of data points in the array.

Our primary finding is that the implementation exhibits faster performance for
lower bitsize variants of the same data type. For instance, int32 demonstrates
a roughly 2x speed improvement compared to int64 for the same number of data
points, while int16 is 2x faster than int32. Remark that these 2x performance
differences are even more pronounced (i.e., a clear 2x) for the parallel
execution. This discrepancy in performance can be attributed to the fact that
reducing the bit-representation by 2x (e.g., int32 vs. int64) allows for a 2x
increase in the number of values that fit in the CPU’s SIMD registers. This
utilization of SIMD registers is an essential part of the ArgMinMax code base
and results in fewer read (memcpy) instructions, which impacts performance
since ArgMinMax is primarily bound by memory access.

Our second key finding emphasizes the benefits of implementing multithreading.
On average, multithreading leads to an impressive 7x performance improvement
on the benchmarking computer. Notably, when extrapolating the linear trend of
the int64 data type, tsdownsample demonstrates the capability to downsample
data at a rate of 45 GB/s (i.e., 8 GB / 0.177 s).

##  5 Conclusion

Time series visualization plays a crucial role in exploratory data analysis,
particularly as datasets continue to grow in size. To enable scalable line-
chart visualization, downsampling has emerged as a well-established technique.
However, we have observed a need for a convenient and high-performance time
series downsampling solution within the Python landscape, facilitating the
integration of downsampling capabilities into widely used Python visualization
packages. To address this need, we introduce tsdownsample, a Python library
that offers a convenient interface to leading precompiled downsampling
algorithms, harnessing highly optimized underlying Rust code. A key aspect of
achieving high performance in tsdownsample involved optimizing argmin and
argmax operations using SIMD instructions and leveraging multithreading. The
runtime feature set detection enables the selection of the most optimal
implementation based on the CPU, allowing for the distribution of a single
binary that can cater to multiple CPU feature sets within the same
architecture. Benchmark results confirm the critical role played by both SIMD
optimizations and multithreading in achieving the impressive performance of
tsdownsample. We firmly believe that tsdownsample represents a significant
advancement in delivering high-performance time series downsampling
capabilities to the Python ecosystem. This advancement is further illustrated
by the adoption of tsdownsample in the plotly-resampler tool, solidifying its
position within the Python community.

## Acknowledgements

The authors thank Martijn Courteaux and Tom Windels for having fruitful
discussions on writing efficient low-level code.

## References

  * [1] W. Aigner, S. Miksch, W. Müller, H. Schumann, C. Tominski, Visualizing time-oriented data—a systematic view, Computers & Graphics 31 (3) (2007) 401–409.
  * [2] B. Shneiderman, The eyes have it: A task by data type taxonomy for information visualizations, in: The craft of information visualization, Elsevier, 2003, pp. 364–371.
  * [3] J. Walker, R. Borgo, M. W. Jones, Timenotes: a study on effective chart visualization and interaction techniques for time-series data, IEEE transactions on visualization and computer graphics 22 (1) (2015) 549–558.
  * [4] N. Bikakis, Big data visualization tools, arXiv preprint arXiv:1801.08336 (2018).
  * [5] P. Godfrey, J. Gryz, P. Lasek, N. Razavi, Interactive visualization of big data, in: Beyond Databases, Architectures and Structures. Advanced Technologies for Data Mining and Knowledge Discovery: 12th International Conference, BDAS 2016, Ustroń, Poland, May 31-June 3, 2016, Proceedings 11, Springer, 2016, pp. 3–22.
  * [6] S. Steinarsson, Downsampling time series for visual representation, Ph.D. thesis (2013).
  * [7] U. Jugel, Z. Jerzak, G. Hackenbroich, V. Markl, M4: a visualization-oriented time series data aggregation, Proceedings of the VLDB Endowment 7 (10) (2014) 797–808.
  * [8] J. Van Der Donckt, J. Van Der Donckt, E. Deprost, S. Van Hoecke, Plotly-resampler: Effective visual analytics for large time series, in: 2022 IEEE Visualization and Visual Analytics (VIS), IEEE, 2022, pp. 21–25.
  * [9] Z. Liu, J. Heer, The effects of interactive latency on exploratory visual analysis, IEEE transactions on visualization and computer graphics 20 (12) (2014) 2122–2131.
  * [10] B. Raskin, N. Aggarwal, The billion data point challenge: Building a query engine for high cardinality time series data (Dec 2018).
URL uber.com/billion-data-point-challenge

  * [11] J. Davi Paganini, Downsampling in the database: How data locality can improve data analysis (Feb 2023).
URL timescale.com/blog/downsampling-in-the-database-how-data-locality-can-
improve-data-analysis/

  * [12] J. Van Der Donckt, J. Van Der Donckt, M. Rademaker, S. Van Hoecke, Minmaxlttb: Leveraging minmax-preselection to scale lttb, arXiv preprint arXiv:2305.00332 (2023).
  * [13] J. Van Der Donckt, J. Van Der Donckt, M. Rademaker, S. Van Hoecke, Data point selection for line chart visualization: Methodological assessment and evidence-based guidelines, arXiv preprint arXiv:2304.00900 (2023).
  * [14] R. Vink, S. de Gooijer, A. Beedie, J. van Zundert, G. Hulselmans, M. E. Gorelli, C. Grinstead, M. Santamaria, D. Heres, ibENPC, J. Leitao, M. van Heerden, C. Jermain, R. Russell, C. Pryer, A. G. Castellanos, J. Goh, M. Wilksch, illumination k, O. Borchert, J. Peek, M. Conradt, L. Brannigan, Y. R. Tan, elbaro, N. Stalder, Robert, S. H. Welling, A. Gregory, pola-rs/polars: Python polars 0.17.3 (Apr. 2023).  doi:10.5281/zenodo.7832012.
URL https://doi.org/10.5281/zenodo.7832012

  * [15] M. Cornea, Ieee 754-2008 decimal floating-point for intel® architecture processors, in: 2009 19th IEEE Symposium on Computer Arithmetic, IEEE, 2009, pp. 225–228.
  * [16] H. Heyman, L. Brandefelt, A comparison of performance & implementation complexity of multithreaded applications in rust, java and c++ (2020).
  * [17] D. Beazley, Understanding the python gil, in: PyCON Python Conference. Atlanta, Georgia, 2010.
  * [18] J. D. Hunter, Matplotlib: A 2d graphics environment, Computing in science & engineering 9 (03) (2007) 90–95.

## Current code version

Nr. |  Code metadata description |  Please fill in this column
---|---|---
C1 |  Current code version |  v0.1.2
C2 |  Permanent link to code/repository used for this code version |  https://github.com/predict-idlab/tsdownsample/releases/tag/v0.1.2
C3 |  Permanent link to Reproducible Capsule |  Not available
C4 |  Legal Code License |  MIT
C5 |  Code versioning system used |  git
C6 |  Software code languages, tools, and services used |  Python, Rust
C7 |  Compilation requirements, operating environments & dependencies |  cargo is used to manage the Rust dependencies and compilation. PyO3/maturin is used as build backend to generate Python bindings for the Rust code.
C8 |  If available Link to developer documentation/manual |  Not available
C9 |  Support email for questions |<EMAIL_ADDRESS>
Table 2: Code metadata (mandatory)

## Current executable software version

Nr. |  (Executable) software metadata description |  Please fill in this column
---|---|---
S1 |  Current software version |  v0.1.2
S2 |  Permanent link to executables of this version |  https://github.com/predict-idlab/tsdownsample
S3 |  Permanent link to Reproducible Capsule |  Not available
S4 |  Legal Software License |  MIT
S5 |  Computing platforms/Operating Systems |  For example Android, BSD, iOS, Linux, OS X, Microsoft Windows, Unix-like , IBM z/OS, distributed/web based etc.
S6 |  Installation requirements & dependencies |  Python 3.7+, Rust nightly
S7 |  If available, link to user manual - if formally published include a reference to the publication in the reference list |  Not available
S8 |  Support email for questions |<EMAIL_ADDRESS>
Table 3: Software metadata (optional)