Olivier Beaumont

Lionel Eyraud-Dubois

Alena Shilova

Efficient Combination of Rematerialization and Offloading for Training DNNs

Rematerialization and offloading are two well known strategies to save memory during the training phase of deep neural networks, allowing data scientists to consider larger models, batch sizes or higher resolution data. Rematerialization trades memory for computation time, whereas Offloading trades memory for data movements. As these two resources are independent, it is appealing to consider the simultaneous combination of both strategies to save even more memory. We precisely model the costs and constraints corresponding to Deep Learning frameworks such as PyTorch or Tensorflow, we propose optimal algorithms to find a valid sequence of memory-constrained operations and finally, we evaluate the performance of proposed algorithms on realistic networks and computation platforms. Our experiments show that the possibility to offload can remove one third of the overhead of rematerialization, and that together they can reduce the memory used for activations by a factor 4 to 6, with an overhead below 20%.

Introduction

A general trend in machine learning, especially in NLP today, is to increase the size of models to improve accuracy. For example, NLP models like BERT [START_REF] Devlin | Bert: Pre-training of deep bidirectional transformers for language understanding[END_REF] or GPT-3 [START_REF] Radford | Language models are unsupervised multitask learners[END_REF] have billions of parameters and vision networks like EfficientNet [START_REF] Tan | Efficientnet: Rethinking model scaling for convolutional neural networks[END_REF], AmoebaNet [START_REF] Asif | Amoebanet: An sdn-enabled network service for big data science[END_REF] or SENET [START_REF] Hu | Squeeze-and-excitation networks[END_REF] routinely have a hundred million parameters, which raises difficult memory issues on architectures like GPUs. Moreover, it was observed [START_REF] Pudipeddi | Training large neural networks with constant memory using a new execution algorithm[END_REF] that a sufficiently large batch-size is required to obtain good and fast convergence. In this paper we focus on the strategies for reducing the memory usage of activations that can be implemented when considering a single computation node, typically a single GPU. These techniques can of course be combined with the different possible parallelizations of the training phase, like data parallelism [START_REF] Zinkevich | Parallelized stochastic gradient descent[END_REF][START_REF] Paine | Gpu asynchronous stochastic gradient descent to speed up neural network training[END_REF][START_REF] Goyal | Accurate, large minibatch sgd: Training imagenet in 1 hour[END_REF], model parallelism [START_REF] Narayanan | PipeDream: generalized pipeline parallelism for DNN training[END_REF][START_REF] Zhan | Pipe-torch: Pipeline-based distributed deep learning in a gpu cluster with heterogeneous networking[END_REF][START_REF] Narayanan | Memory-efficient pipeline-parallel dnn training[END_REF], filter or kernel [START_REF] Nikoli Dryden | Channel and filter parallelism for large-scale cnn training[END_REF] or image parallelisms [START_REF] Nikoli Dryden | Improving strong-scaling of cnn training by exploiting finer-grained parallelism[END_REF]. The combination with parallel strategies is however out of the scope of the present paper and our goal here is to model, analyze the complexity and understand how to benefit simultaneously from rematerialization and offloading.

Rematerialization [START_REF] Gruslys | Memory-efficient backpropagation through time[END_REF][START_REF] Chen | Training deep nets with sublinear memory cost[END_REF][START_REF] Kukreja | Backpropagation for long sequences: beyond memory constraints with constant overheads[END_REF][START_REF] Beaumont | Optimal checkpointing for heterogeneous chains: how to train deep neural networks with limited memory[END_REF] is a technique that consists in deleting from memory some of the activations computed during the forward phase as soon as they are no longer needed for the forward phase. As these activations are required again during the backward phase, they are then recomputed from other activations kept in memory. The idea is therefore to trade additional computations for lower memory requirements. The objective, in the context of rematerialization, is to compute a sequence of valid operations (i.e. such that the inputs of any operation are actually in memory at the time of the operation) consisting of elementary forward, backward and delete operations, and which has a minimum execution time among all sequences that satisfy a given memory constraint. Offloading [START_REF] Rhu | vdnn: Virtualized deep neural networks for scalable, memory-efficient neural network design[END_REF][START_REF] Shriram | Dynamic memory management for gpu-based training of deep neural networks[END_REF][START_REF] Tung D Le | Tflms: Large model support in tensorflow by graph rewriting[END_REF][START_REF] Beaumont | Optimal gpu-cpu offloading strategies for deep neural network training[END_REF] is a technique that has the same objective, i.e. to produce a sequence of operations satisfying a given memory constraint. However, in the case of offloading, each computation operation is executed exactly once.The idea here is therefore to trade transfers on the PCI bus between the CPU and the GPU for lower memory requirements. The objective, in the context of offloading, is therefore to compute a sequence of valid operations consisting of elementary forward, backward, offload and prefetch operations, and which has a minimum execution time among all sequences that fulfill a given memory constraint.

While rematerialization and offloading share the common goal of limiting memory consumption, they rely on different techniques and achieve this result by consuming independent resources, using either more computation or more transfers between a GPU and a CPU. In this paper we solve the problem of finding the optimal sequence combining rematerialization and offloading, which has not been addressed in the literature to the best of our knowledge. We discuss related work in Section 2 and model and notations in Section 3. The algorithmic solution to combine optimally checkpointing and offloading is presented in Section 4 and experiments to assess the advantage of using both offloading and rematerialization, based on a large range of neural networks and architecture parameters, are presented in Section 5.

Related work

Parallelism-based memory optimization A first approach to control memory consumption is to use parallelization. Two main sources of memory consumption in the training phase for deep neural networks are, on the one hand, the storage of the network weights and, on the other hand, the storage of activations. To limit the storage cost of activations (for a fixed batch size), data parallelism [START_REF] Zinkevich | Parallelized stochastic gradient descent[END_REF][START_REF] Paine | Gpu asynchronous stochastic gradient descent to speed up neural network training[END_REF][START_REF] Goyal | Accurate, large minibatch sgd: Training imagenet in 1 hour[END_REF] is the method of choice. Each batch is divided into mini-batches which are in turn processed in parallel on different computation resources. In this scheme, however, weights are replicated on each resource and collective communications are required to synchronize weights between different resources. When models become large, these collective communications can hinder the scalability of the data parallel approach.

On the contrary, model parallelism [START_REF] Dean | Large scale distributed deep networks[END_REF][START_REF] Huang | Gpipe: Efficient training of giant neural networks using pipeline parallelism[END_REF][START_REF] Narayanan | PipeDream: generalized pipeline parallelism for DNN training[END_REF][START_REF] Zhan | Pipe-torch: Pipeline-based distributed deep learning in a gpu cluster with heterogeneous networking[END_REF] is used to distribute weights over different available resources. Weights of a neural network are then distributed on different computation resources, and communications happen between consecutive layers assigned to different resources. Nevertheless, this approach also requires sophisticated strategies [START_REF] Narayanan | Memory-efficient pipeline-parallel dnn training[END_REF] to save fewer activations and models in memory. The data can also be split along other dimensions. In [START_REF] Nikoli Dryden | Improving strong-scaling of cnn training by exploiting finer-grained parallelism[END_REF], large images are split into smaller images, the network is trained on these small images (augmented by a halo) and extra communications are needed in order to synchronize parameter updates. A similar strategy was proposed to use channel and filter parallelism in [START_REF] Nikoli Dryden | Channel and filter parallelism for large-scale cnn training[END_REF]. All these parallel strategies can be combined with the single node memory optimization strategies considered in the present paper.

Rematerialization-based strategies Rematerialization techniques consist in keeping only a subset of activations in memory during the whole training phase, while the others are dynamically recomputed at runtime when they are needed for the backpropagation. This allows to explore a tradeoff between memory usage and computational cost. Rematerialization was first considered in the context of Automatic Differentiation (AD) [START_REF] Griewank | On automatic differentiation[END_REF] under the name "checkpointing". In the context of AD, computational networks consist in very long and homogeneous chains, in which the forward activation corresponding to the i-th stage of the chain has to be kept into memory until the i-th backward stage. Rematerialization algorithms consist in determining in advance which forward checkpoints should be kept into memory and which ones should be recomputed from stored checkpoints when performing the backward phase. In the case of homogeneous chains closed form formulas providing the exact position of checkpoints were proposed in [START_REF] Griewank | Algorithm 799: Revolve: an implementation of checkpointing for the reverse or adjoint mode of computational differentiation[END_REF].

The use of rematerialization strategies was advocated for DNN in several papers [START_REF] Gruslys | Memory-efficient backpropagation through time[END_REF][START_REF] Chen | Training deep nets with sublinear memory cost[END_REF][START_REF] Kukreja | Backpropagation for long sequences: beyond memory constraints with constant overheads[END_REF]]. An implementation of a simple periodic strategy [START_REF] Chen | Training deep nets with sublinear memory cost[END_REF] was provided in PyTorch [START_REF] Paszke | Automatic differentiation in pytorch[END_REF] 1 , but it is only optimal for the restricted case of homogeneous chains and when activation discards are done only before the loss. However, DNN models are in general more complicated. More recently, this strategy has been extended to general heterogenous networks represented as sequential networks in [START_REF] Beaumont | Optimal checkpointing for heterogeneous chains: how to train deep neural networks with limited memory[END_REF] 2 .

The case of non sequential networks was also addressed in the literature [START_REF] Kirisame | Dynamic tensor rematerialization[END_REF][START_REF] Kumar | Efficient rematerialization for deep networks[END_REF][START_REF] Jain | Checkmate: Breaking the memory wall with optimal tensor rematerialization[END_REF][START_REF] Kusumoto | A graph theoretic framework of recomputation algorithms for memory-efficient backpropagation[END_REF][START_REF] Feng | Optimal gradient checkpoint search for arbitrary computation graphs[END_REF], mostly under the simplifying assumption that discarding computed activations is only possible before the computation of the loss function [START_REF] Kusumoto | A graph theoretic framework of recomputation algorithms for memory-efficient backpropagation[END_REF][START_REF] Feng | Optimal gradient checkpoint search for arbitrary computation graphs[END_REF]. In [START_REF] Kirisame | Dynamic tensor rematerialization[END_REF] the authors proposed an algorithm dealing with dynamic data flows, though based on an heuristical approach. The authors of [START_REF] Jain | Checkmate: Breaking the memory wall with optimal tensor rematerialization[END_REF] considered minimization of the runtime with the use of Integer Linear Programming. Even though this approach is suitable for arbitrary computational graphs with heterogeneous costs, solving this ILP is computationally very expensive even for intermediate depth networks.

Offloading-based strategies Offloading [START_REF] Rhu | vdnn: Virtualized deep neural networks for scalable, memory-efficient neural network design[END_REF][START_REF] Shriram | Dynamic memory management for gpu-based training of deep neural networks[END_REF] is an alternative approach which consists in offloading some of the forward activations from the memory of the GPU to the memory of the CPU, which is expected to be much larger. In [START_REF] Rhu | vdnn: Virtualized deep neural networks for scalable, memory-efficient neural network design[END_REF], the authors proposed a simple mechanism of memory virtualization, that introduces unnecessary synchronizations between data transfers and computations of later forward activations. This technique was improved in [START_REF] Shriram | Dynamic memory management for gpu-based training of deep neural networks[END_REF]. In both papers, the algorithmic strategies to decide which activations to offload into the main memory are relatively straightforward, and basically consist in sending to the CPU either all activations or only those of convolutional layers. In [START_REF] Beaumont | Optimal gpu-cpu offloading strategies for deep neural network training[END_REF], the complexity of the offloading problem was analyzed and several NP-completeness results were provided, depending whether activations can be offloaded partially or not.

In order to reduce the overhead incurred by the communications, the use of compression was advocated in [START_REF] Rhu | Compressing dma engine: Leveraging activation sparsity for training deep neural networks[END_REF]. Memory virtualization was further considered in [START_REF] Meng | Training deeper models by gpu memory optimization on tensorflow[END_REF][START_REF] Tung D Le | Tflms: Large model support in tensorflow by graph rewriting[END_REF][START_REF] Huang | Swapadvisor: Pushing deep learning beyond the gpu memory limit via smart swapping[END_REF][START_REF] Zhang | Efficient memory management for gpu-based deep learning systems[END_REF]. Finally, the speed-centric rematerialization from [START_REF] Chen | Training deep nets with sublinear memory cost[END_REF] enhanced with memory-centric rematerialization (discards activations of every segment all the time) was combined with the simple offloading approach from [START_REF] Rhu | vdnn: Virtualized deep neural networks for scalable, memory-efficient neural network design[END_REF] in [START_REF] Wang | Superneurons: Dynamic gpu memory management for training deep neural networks[END_REF]. In addition, the authors in [START_REF] Rajbhandari | Zero: Memory optimizations toward training trillion parameter models[END_REF] also used rematerialization of [START_REF] Chen | Training deep nets with sublinear memory cost[END_REF] with a possibility of further offloading saved checkpoints to the CPU if rematerialization only is not enough to perform training under memory constraint.

Another approach that combines recomputations and data offloading from GPU memory to CPU memory was proposed in [START_REF] Pudipeddi | Training large neural networks with constant memory using a new execution algorithm[END_REF]. This approach is especially useful in the case where the size of the activations is small compared to the size of the model, which is the case in some NLP models. In this case, the network weights are offloaded to the CPU memory, that serves as a parameter-server host. In this paper, we focus on the case where the model can reside in the memory of a single GPU and the activations should be moved to the CPU memory. The analysis of the combination of these two approaches is out of the scope of this paper.

Model and notations

Forward, backward and offloading elementary operations

In what follows, our goal is to process the chain described by Figure 1 that corresponds to the processing of a mini-batch. This chain depicted by Figure 1 is defined by the set of operations F and B . The default forward operations F requires only the activation produced by the previous layer F -1 (or the input mini-batch in the case of F 1 ), whereas the default backward operation B require the preceding gradient δ , the activation a -1 and the intermediate generated data ā necessary for the backward step. The distinction between a and ā is crucial in practice, since it allows to consider computation graphs which have a sequential structure, but are not purely sequential. In this setting, it is indeed possible for F to represent a complex operation (any Direct Acyclic Graph of layers), and thus ā contains all the intermediate activations produced by F and necessary to compute B , whereas a only contains the output activation of F , that will be used by F +1 . In general, a can therefore be much smaller than ā .

In order to manipulate different types of activations, we define in Table 1 a set of generic forward and backward operations, that can be directly translated in automatic differentiation packages such as tf.GradientTape for Tensorflow or torch.autograd for PyTorch. Operations are defined both by the type of operands they use and produce (a or ā ) and by the treatment they make of the input data after processing (deleted or kept in memory). Each operation is associated to a processing time and a memory overhead that takes place during the computation. Both the processing time and the memory overhead only depend on whether the operation is a forward or a backward task.

F 1 F 2 • • • F L-1 F L loss B 1 B 2 B 3 • • • B L B loss a 0 a 1 a 2 a L-2 a L-1 a L a L+1 δ L+1 = 1 δ L δ L-1 δ 3 δ 2 δ 1 δ 0 a 0 a 1 a 2 a L-1 a L ā1 ā2 ā3 āL-1 āL āL+1
-1 } {a -1 , ā } u F o f {ā -1 } {ā -1 , ā } F ck Forward and materialize input {a -1 } {a -1 , a } u F o f {ā -1 } {ā -1 , a } F ∅ Forward without saving {a -1 } {a } u F o f {ā -1 } {a } B l Backward step {δ , ā , a -1 } {δ -1 } u B o b {δ , ā , ā -1 } {δ -1 , ā -1 } Table 1: Generic operations available in DL frameworks.
In addition to the operations described in Table 1, we also consider two operations O and P that correspond respectively to offloading to and prefetching from CPU memory. O and P can operate both on a -1 and ā -1 , depending on which type of activation has been produced. Thanks to CUDA streams, it is possible to asynchronously offload and prefetch activations from and to the CPU memory, while independent computations are performed.

Assumptions and limitations

In what follows, our goal is to find a sequence of F all , B , F ck , F ∅ , O and P operations to compute δ 0 from a 0 . Given a bandwidth β between the GPU and the CPU, we can compute the peak memory usage and the total duration of a particular sequence. Our goal is to find a sequence (i) whose memory peak is smaller than the available memory on the GPU M GPU , and (2) whose overall duration is minimal. In order to make the problem tractable, we introduce additional assumptions 1. Offloading before loss: all offloading operations must complete before the computation of the loss and prefetching operations (each prefetching operation being performed only once) start only after the computation of the loss. This is a reasonable assumption since offloading (resp. prefetching) introduce delays and therefore should be performed as early (resp. as late) as possible.

2. Memory persistency holds: if for some layer i activation āi is computed, then no operations related to layers i with i < i take place until B i . This is a common assumption in the rematerialization literature, and it was shown in [START_REF] Beaumont | Optimal checkpointing for heterogeneous chains: how to train deep neural networks with limited memory[END_REF] that in all practical cases, optimal solutions satisfy this assumption (despite a theoretical result showing that there are cases where this is not true)

3. Asynchronous and continuous offloading: Activations should be offloaded/prefetched entirely, but memory allocation and release can be performed in several steps during the communication. This assumption proved to be efficient in [START_REF] Beaumont | Optimal gpu-cpu offloading strategies for deep neural network training[END_REF] to provide a tractable context without degrading the solution quality. As the practical implementation of this multi-step offloading is not trivial, for the experimental evaluation in Section 5 we transform the obtained solutions into solutions in which each communication is done at once.

In Section 4, our goal is to find the optimal sequence of operations under above assumptions. Removing some of these assumptions may therefore improve the solution (by seeking the solution in a larger space), but a priori at the cost of more expensive computations.

1. A first limitation comes from the ability to offload ā activation, which is not easy to implement in frameworks like PyTorch, because it requires modifying and accessing internal structures. Solutions that currently implement offloading, such as VDNN [START_REF] Rhu | vdnn: Virtualized deep neural networks for scalable, memory-efficient neural network design[END_REF][START_REF] Shriram | Dynamic memory management for gpu-based training of deep neural networks[END_REF], therefore require a thorough rewriting of the transfer operations.

2. A second limitation of our work is that it relies on the possibility of linearizing networks, i.e. representing them as a sequence of operations, in which each operation can in turn be arbitrarily complex (and associated to any type of dependences). In practice, most networks are easy to sequentialize. In the case of rematerialization only, optimal solutions based on integer linear programming such as Checkmate [START_REF] Jain | Checkmate: Breaking the memory wall with optimal tensor rematerialization[END_REF] were proposed. However, the cost of integer linear programming limits their use to very small networks (smaller than those considered in Section 5) and their relaxation come without any guarantee. The extension of the proposed approach to complex dependency structures is nevertheless not straightforward.

Combination of offloading and rematerialization

We propose a dynamic programming to optimally solve the problem of finding a sequence containing rematerialization and offloading operations of minimal duration. This algorithm is denoted by pofo, for "Persistent with Offloading during Forward Only". To save space, we will not detail in the main paper all the equations of the dynamic program, which involves a large number of cases. We will focus in the main part of the paper on the intuitions and the general working principle of the dynamic program and refer the reader to Appendix B for detailed derivations and proofs.

Rationale of the different operations

The pofo algorithm is based on a sequence of choices, that consist in deciding, for each 1 ≤ i ≤ L (i) which type of operation F i we are going to use and whether we are going to compute a i or āi (ii) whether we are going to keep the input value a i-1 or āi-1 in the memory of the GPU, offload it in the memory of the CPU or completely delete it from the memory and (iii) how to compute B i . Concerning (i), the size of āi is in general larger than a i . Both can be used to compute F i+1 and B i+1 , but if a i only is kept in memory (either CPU or GPU), then F i will have to be recomputed during the backward phase before B i . Concerning (ii), if we delete a i-1 or āi-1 from memory, then it will be necessary to recompute F i-1 before processing B i , but it will save memory for the whole subsequent sequence F i+1 . . . F L B L . . . B i+1 . If a i-1 or āi-1 is offloaded on CPU memory (which may take some time), the memory will be available for the subsequent sequence, and a i-1 or āi-1 will have to be prefetched before computing B i . The interest of offloading obviously depends on the bandwidth β of the PCI bus. Concerning (iii), if the input activation (a i-1 or āi-1 ) has been deleted from memory, we will have to recompute it, starting from the last kept activation (either a k or āk for some k < i -1). Several rematerialization sequences are possible to do this operation that have potentially different durations and that offer different prefetching possibilities, depending on their memory profile, and the choice of the optimal sequence will depend on the memory state before computing B i , as detailed in Appendix B.

Intuition of the overall scheme and state variables

The dynamic programming solution is built in several steps. We consider separately the forward propagation within which the offloading of activations can be done. The forward propagation is followed then with the loss calculation synchronized with the end of offloading operations. After loss computation, the backward propagation interleaved with the prefetching starts.

The general principle of pofo is the following. During the forward phase (before the computation of the loss), we consider the layers in an increasing order (from 1 to L). At each step, we have to decide which operation to implement among F all , F ck , F ∅ and in the case where the input activation is kept, we have to decide if we add O . Inspired by the ideas of Automatic Differentiation [START_REF] Griewank | Algorithm 799: Revolve: an implementation of checkpointing for the reverse or adjoint mode of computational differentiation[END_REF] (but proposed in the context of a fully homogeneous chain without offloading), our dynamic programming relies on the following remark: any solution can be decomposed into parts, where each part computes layers i to j (j ≥ i), and only the input to layer i is stored in memory. Starting from a layer i, whose input needs to be saved in memory, pofo needs to decide whether the input will be saved with F i ck or F i all . It also needs to decide the index j of the end of the corresponding part in the solution, so that j + 1 starts the next part of the sequence, and its input will be the first activation kept in memory after i. A recursive call to the sub-chain starting from j allows to obtain the corresponding running time. Between the two saved activations of layers i and j, a sequence of F k ∅ takes place. As forwards F i ck and F i all correspond to two different behaviors, they generate two different cases of the dynamic programming that are considered in Appendix B.1. In addition, an offloading decision needs to be made for the input of layer i, either store it in the memory of the GPU or offload it to the CPU; this decision impacts the memory available for the rest of the chain and idle times due to communications.

In order to use dynamic programming, we need to define a set of variables that describes the state of the system at any instant. This set of variables should be as small as possible, since it has a direct influence on the size of the data structure and on the computing time to solve the dynamic program. On the other hand, these variables must be chosen wisely and they must contain enough information to make decisions for subsequent layers and be updated according to these decisions. To evaluate correctly memory constraints and to compute idle times, it is important to know how memory may vary between layers i and j with saved inputs: from some minimal memory occupation till maximal memory occupation. Both can be described for the forward and the backward phases using only three state variables: in addition of the index i of the last layer with saved input discussed above, we use A i , ∆ Fi and ∆ Bi , in addition to a boolean variable x.

• A i denotes the total GPU memory occupied by the saved values among a 0 , . . . , a i-2 and ā1 , . . . , āi-2 that are not transferred to the CPU.

• ∆ Fi denotes the amount of data from a 0 , . . . , a i-2 and ā1 , . . . , āi-2 that the schedule still needs to offload after F i-1 .

• ∆ Bi denotes the amount of data from a 0 , a 1 , . . . , a i-2 and ā1 , . . . , āi-2 that the schedule should prefetch before starting B i-1 .

• x is a boolean specifying whether the input of layer i is

a i-1 (if x = 0) or āi-1 (if x = 1).
We denote as

I x i = (1 -x)a i-1 + xā i-1
the size of the input activation for layer i. With these state variables, we can compute M x Fi , the memory on the GPU after executing F i-1 , and M x Bi , the memory on the GPU before executing B i-1 (excluding δ i-1 ) are given by M x Fi = A i + ∆ Fi + I x i and M x Bi = A i + ∆ Bi + I x i . These memory values represent the maximal memory occupation between layers i and j for the forward and backward phases respectively, and are used for computing idle times when overlapping communications with computations, according to Lemma 1. On the other hand, to estimate the feasibility of the scheduled operations, we also need to know the memory after everything has been offloaded, which is either A i + I x i (when I x i stays on the GPU) or A i (when I x i is sent to the CPU), which we use to obtain the maximal memory available on the GPU.

The updates of the state variables follow simple rules described below

• if no new data is offloaded (prefetched), then ∆ Fi (∆ Bi ) is constantly decreasing from index i to index j at speed β, until reaching zero;

• if new data has to be offloaded (prefetched), then the data size is added into ∆ Fi (∆ Bi );

• A i is updated if new data is saved on the GPU (without being later offloaded to the CPU).

Our goal is to find the sequence of operations that minimizes the overall execution time. However, both rematerialization and offloading can induce extra time with respect to the execution of the sequence F 1 . . . F L B L . . . B 1 which can be computed given infinite memory on the GPU:

• During the forward phase, offloading helps to keep GPU memory low, but, if the transfers to the CPU are not fast enough due to limited bandwidth, some idle time might occur, waiting for enough memory to be freed by offloading. This is analyzed in details in Appendix B.1.

• At the interface between the forward and the backward phase, we need to enforce that offloading and prefetching are well synchronized by the computation of the loss, what can in turn introduce some idle time on the GPU. This will be analyzed in Appendix B.2.

• During the backward phase (after the computation of the loss), there might be two sources of delays: recomputations of discarded activations and idle times induced by prefetching. Indeed, some prefetching operations may not be completed by the time the activation is needed, if they are delayed because of memory constraints. In fact, these two sources of extra times are not independent, since a longer rematerialization sequence with lower memory needs might allow more overlapping of prefetching with computations and thus avoiding idle time later. We propose an auxiliary dynamic program that includes the computation of intermediate idle times in each recursive call to find the best schedule under prefetching. The idle times for prefetching are found using Lemma 1. This is analyzed in Appendix B.3.

Finally, we present a lemma that describes the general scheme for optimally overlapping the communications with scheduled operations. Lemma 1 is valid for both prefetching and offloading, as their behavior is symmetrical. Offloading takes place at the beginning of the sequence, making available memory increase at the speed β from its initial value m min . On the contrary, prefetching is done at the end of the sequence, making available memory decrease at the speed β until reaching m min . Performing offloading as soon as possible and prefetching as late as possible allows to have more available memory for the middle of the execution. Note that m min could be equivalently replaced with M GPU -M max , where M max denotes the maximal memory occupation (either M x Fi or M x Bi ).

Lemma 1 Let us consider a fixed sequence of operations, for which the available memory increases (resp. decreases) during execution because of data unloading (resp. prefetching), with its minimum at m min . Let us denote by M S o the memory required to process operation o ∈ S, and d o the distance between o and the end (beginning) of sequence S,i.e. the cumulative duration of operations taking place before (after) operation o. Then, the execution of S needs to be delayed by some idle time

= max max o∈S (M S o -βd o ) -m min β , 0 .
The above mentioned ingredients can be used to build the dynamic programming algorithm pofo, solving optimally the problem in the pseudo-polynomial time. Theorem 1 states this result, while the complete description of the dynamic program and the proof of the theorem can be found in Appendix B. The algorithm is pseudo-polynomial in the memory limit and it relies on discretization of the memory values. In the experiments of Section 5, we systematically use 50 discretization steps (scaling so that M GPU = 50), which shows no discretization effect (we tried more steps to confirm this) and a reasonable complexity, with an execution time below 4 minutes in the worst cases. Since optimization is performed only once before the training starts, we argue that this is fully acceptable.

Theorem 1 Under assumptions of Section 3.2, the problem of finding the minimal processing time for the chain from Figure 1, using operations from Table 1 together with offloading O and prefetching P , under memory limit M GPU discretized with N GPU values and bandwidth β, can be solved optimally with a dynamic programming algorithm with a complexity of O(L 2 N 3 GPU + L 3 N 2 GPU ).

Experiments

Additional heuristics

We have implemented two other sophisticated heuristics for comparison purposes. These heuristics are not bound with the assumptions of Section 3.2, but come without any optimality for the produced sequences.

In the first heuristic, called opportunist, the objective is to use the communication medium as much as possible. We compute the first layer with F all and offload its input and output, ensuring that the memory will remain fully available for the rest of the computation. The next layers are computed with F ∅ until the end of communications. We then start a new communication by performing next layer with F all and offloading its input and output, and so on until the end of the sequence can be entirely performed in memory. This process thus builds groups of layers between two F all operations. We then compute the backward phase for each group using the rematerialization algorithm, and concatenate them (with the necessary prefetches) to obtain the final sequence. Note that opportunist is conceptually close to the implementation from Superneurons [START_REF] Wang | Superneurons: Dynamic gpu memory management for training deep neural networks[END_REF].

The second heuristic is called autocapper. It relies on an internal capper algorithm, that uses an increased memory limit M > M GPU as an additional input. capper computes a pure rematerialization sequence with the limit M , finds the peak memory usage in the sequence and offloads the lowest indexed activation present in GPU memory at that time. This process is repeated until the sequence fits in the memory limit M GPU . autocapper calls capper with 40 values of M , evenly spaced between the target M GPU and M high , the memory required without rematerialization or offloading. For each value M , the resulting sequence from capper is simulated and the best one is kept. Note that, when autocapper does not perform recomputations, it behaves as the Greedy heuristic in [START_REF] Beaumont | Optimal gpu-cpu offloading strategies for deep neural network training[END_REF].

Experimental setting

We measured running times and memory occupation of several networks from PyTorch torchvision package: resnet, densenet and inception, with a batch size of 16 and images of 500×500 pixels. We also experimented with other batch and image sizes, and obtained very similar results. Time measurements were performed on a NVidia Tesla V100 GPU. We also measured the bandwidth obtained when transferring PyTorch tensors from and to the GPU and obtained 12GB/s. The simulation results presented here were obtained using 4 cores of a 24-core Haswell Intel R Xeon R E5-2680 v3 at 2,5 GHz, with 128GB of memory, and used about one hour of total computation. We have modified the open-source rematerialization framework rotor, developed by Inria, to handle offloading, and we have implemented our algorithms in this framework.

We consider 5 algorithms in total: our dynamic program pofo, the optimal rematerialization-only algorithm from [START_REF] Beaumont | Optimal checkpointing for heterogeneous chains: how to train deep neural networks with limited memory[END_REF], the optimal offload-only approach (DynProg heuristic from [START_REF] Beaumont | Optimal gpu-cpu offloading strategies for deep neural network training[END_REF]), and both heuristics autocapper and opportunist. We use the sequential time as a reference: it is the time that it would take to process the forward and backward phases with infinite memory, equal to the sum of forward and backward times of all operations. For each network, we compute the highest and lowest memory requirements (denoted respectively M high and M low ): M high is obtained with the sequential approach, while M low is obtained by recomputing everything from the beginning at each step of the backward phase. We can thus explore the whole range of achievable memory sizes for this network, by considering values within the interval [M low , M high ]: for a given ratio α ∈ [0, 1], the memory limit M GPU is set to (1 -α)M low + αM high . For the algorithms which perform offloading, we also vary the bandwidth value β, from 12GB/s (corresponding to a realistic scenario) to 36GB/s (corresponding to possible future improvements in communication capabilities). The results are shown on two plots: Figure 2 presents the behavior of the algorithms for fixed bandwidth and varying memory constraint, while Figure 3 explores the effects of increasing the bandwidth for a given memory limit. We draw the following conclusions:

• When memory is large (α = 0.75), offloading is more effective than rematerialization: with β = 12GB/s, there is time to offload and prefetch enough data to run the whole chain. However, for smaller memory limits, the offload-only policy exhibits much worse performance than rematerialization (unless more bandwidth is available, see Figure 3).

• Unless in some cases with small memory, the opportunist policy achieves slightly better performance than pure rematerialization, but is sometimes unable to produce a solution when the memory limit is too low. On the other hand, the more sophisticated autocapper algorithm obtains very good performance for memory ratios above 0.5.

• The pofo algorithm successfully combines the advantages of both rematerialization and offloading, and consistently outperforms both of them. When the memory limit is high enough, the optimization problem is relatively easy, and autocapper and pofo achieve similar performance. In some cases autocapper is marginally better than pofo, as pofo produces solutions under assumptions of Section 3.2, while partial memory releases (assumption (iii)) are difficult to implement in practice. For lower memory limits, the more optimized pofo algorithm is able to produce much better sequences.

• As can be seen on Figure 3, the offload-only approach works perfectly when the bandwidth is high. Indeed, if there is enough communication capability to offload all the data while it is produced, it is possible to avoid recomputations without inducing idle time. Compared to pure offloading techniques, our pofo solution allows to train models with large activations with cheaper communication links.

Conclusion

In this paper, we formalize the problem of the optimal combination of rematerialization and offloading, which are two classical strategies for coping with memory limitations on a GPU. We show that the optimal solution can be computed using dynamic programming, in a few seconds or a few minutes for very deep networks.

From experiments, we show that the combination of offloading and rematerialization is very efficient and allows, in many cases, to transparently perform training with 4 to 6 times less memory, at the cost of a 10-20% time overhead. Among the perspectives of this work, one could be interested in minimizing the energy consumption under both memory and throughput constraints. This study would enable to quantify the relative effects of data transfers and redundant computations from a different viewpoint. Another perspective of this work would be to adapt the dynamic programming approach to the case where we consider the possibility of offloading not only the activations, but also the weights of the network, which is necessary in the case of recent very large NLP networks in particular and which is used, for example, on the basis of heuristic optimization in Deepspeed3 . The dynamic programming approach also seems to be adapted to this new context.

A Proof of Lemma 1

Let us concentrate on prefetching, while the case with offloading is done analogously. The idle time is not zero when it is not possible to overlap entirely communications and computations because of the memory.

During prefetching the available memory (excluding the memory needed for operations in S) is decreasing at speed In this section, we detail the dynamic programming equations whose intuition has been given in Section 4.

β
In order to proceed, one should know the minimal memory requirements for each operation. Given the chain between i and j, taking into account the data dependencies from Table 1 and the fact that the global input of the chain

I x i = (1 -x)a i-1 + xā i-1
(it is equal to a i-1 when x = 0 and āi-1 when x = 1) is already stored in the memory (and counted separately in all equations) and the current gradient value should be stored at any time, then

• M i,j F i ∅ = δ j + a i + o f i ; • ∀k = i : M i,j F k ∅ = δ j + a k-1 + a k + o f k ; • M i,j F i all = δ j + āi + o f i ; • ∀k = i : M i,j F k all = δ j + a k-1 + āk + o f k ; • M i,j Bi = āi + δ i + δ i-1 + o b i ; • ∀k = i : M i,j B k = a k-1 + āk + δ k + δ k-1 + o b k .

B.1 Forward phase

In this section, we provide the analysis of the forward phase,i.e. for all operations that take place before the computation of the loss. We consider the subchain that starts from layer i and assume that its input I x i is saved by F i ck or F i all . To denote the total duration of the execution from F i to B i , we use OC x (i, A i , ∆ Fi , ∆ Bi ). We also distinguish between OC F all x (resp. OC F ck x ) that represents the duration of the chain execution when the first operation is constrained to be F i all (resp. F i ck ). By construction, the first operation cannot be F i ∅ since we assume that the input of the chain must not be discarded. Case 1: F i all is the first operation If the first operation is F i all , then in the optimal schedule the next operation cannot be F i+1 ∅ . Indeed, as F i all stores āi together with the input I x i , then this āi should stay in the memory until the backward operation B i . At the same time, as F i+1 ∅ discards its input after its execution, then each time when we need to recompute a i we need to re-execute F i . Therefore, if the next operation is F i+1 ∅ , the overall duration does not change if the first instance of F i all is replaced by F i ck , while F i all replaces the last F i that computes a i in the sequence. Simultaneously, this transformed sequence keeps less data in memory, which contradicts the optimality of F i all being followed by F i+1 ∅ .

offload input Since we assume memory persistency,i.e. that āi stays in the memory until the backward operation B i , the sequence after F i all and before B i can be computed with a recursive call to the dynamic programming. After this sequence, B i can be directly executed, having all the necessary input already stored in device memory. This shows that OC F all x (i, A i , ∆ Fi , ∆ Bi ) can be computed as

no offload v = 1 v = 2 new A v A i A i + I x i new ∆ v F max{∆ Fi + I x i -D F , 0} max{∆ Fi -D F , 0} new ∆ v B max{∆ Bi -D B , 0} + I x i max{∆ Bi -D B , 0}
OC F all x (i, A i , ∆ Fi , ∆ Bi ) = u Fi + u Bi + min v=1,2 OC x=1 (i + 1, A v i+1 , ∆ v Fi+1 , ∆ v Bi+1 ) + F + B ,
where the values F and B represent the idle time required for communications in the forward and backward respectively, and can be determined using Lemma 1. Taking into account that the maximal memory occupation for F i all and B i are given by M x Fi and M x Bi respectively, we obtain

F = max   M i,L F i all -M GPU + M x Fi β , 0   and B = max M i,L Bi -M GPU + M x Bi β , 0 .
We also need to compute the new state variables

A v i+1 , ∆ v Fi+1 , ∆ v Bi+1
, where the value of v indicates whether the input I x i is offloaded or not. During the forward step, the communication link is able to offload a quantity of data given by D F = (u Fi + F )β; similarly, the data that can be prefetched during the backward step is given by D B = (u Bi + B )β. These values can be used to obtain the new state variables, as described in Table 2.

The final sequence is valid if memory limit is not violated. So if

M i,L F i all + I x i > M GPU -A i or M i,L Bi + I x i > M GPU -A i then we set OC F all x (i, A i , ∆ Fi , ∆ Bi ) = ∞. Case 2: F i
ck is the first operation After an F i ck operation, any forward operation or offload operation is possible. Let us assume that the next saved activation is a j for some j ≥ i, which implies that after F i ck there is a sequence of F k ∅ for i < k ≤ j. Due to memory persistency, it also implies that the sequence processing layers from j + 1 till L can be obtained recursively, looking at OC x=0 (j + 1, A v j+1 , ∆ v Fj+1 , ∆ v Bj+1 ) (once again, v denotes whether the input I x i is offloaded or not). Afterwards, processing B k for i < k ≤ j requires recomputing forward operations from a i-1 (or āi-1 ). Since no offloading operation can be performed after the computation of the loss, this corresponds to computing a rematerialization sequence between layers i and j that should have ∆ Bi prefetched by the end. C ∆ (i, j, A i + I x i , ∆ Bi ) is used to compute the optimal duration to process layers from i to j, having A i + I x i as a cumulative storage of all activations at the GPU, and ∆ Bi that corresponds to the data that has to be prefetched. We show in Section B.3 how to compute C ∆ (i, j, A, ∆). This yields the formula for

OC F ck x (i, A i , ∆ Fi , ∆ Bi ) OC F ck x (i, A i , ∆ Fi , ∆ Bi ) = min i≤j≤L-1 ( j k=i u F k + min v=1,2 OC x=0 (j + 1, A v j+1 , ∆ v Fj+1 , ∆ v Bj+1 ) +C ∆ (i, j, A i + I x i , ∆ Bi ) + F
), where F is computed with the help of Lemma 1

F = max   max k=i,...,j (M i,L F k ∅ -β k-1 h=i u F h ) -M GPU + M x Fi β , 0   .
Similarly to the Case 1, the values of the new state variables A v , ∆ v F and ∆ v B can be found using the formulas in Table 2. The only difference is the possible amount of offloaded data and prefetched data: in this case

D F = ( j k=i u F k + F )β and D B = C ∆ (i, j, A, ∆ Bi )β.
The final sequence is valid if memory limit is not violated. Thus, if for some k ≥ i we have M

i,L F k ∅ + I x i > M GPU -A i when I x i is not offloaded (or when k = i) or M i,L F k ∅ > M GPU -A i when I x i is offloaded, then we set OC F ck x (i, A i , ∆ Fi , ∆ Bi ) = ∞.
Combining everything together Therefore, OC x (i, A i , ∆ Fi , ∆ Bi ) can be computed as

OC x (i, A i , ∆ Fi , ∆ Bi ) = min OC F all x (i, A i , ∆ Fi , ∆ Bi ) OC F ck x (i, A i , ∆ Fi , ∆ Bi ) (1)
B.2 Loss: how to concatenate forward and backward phases

In the previous section, we have shown how to compute the optimal duration of the sequence, using dynamic programming with OC x (i, A i , ∆ Fi , ∆ Bi ). This dynamic program finds the solution through recursive calls to smaller sub-chains, until reaching the subchain consisting of only loss computation. We further represent the loss computation with F L+1 and B L+1 operations. The loss should be computed with F L+1 all and B L+1 , so that this case is similar to OC F all x and

OC x (L + 1, A i , ∆ Fi , ∆ Bi ) = u F L+1 + u B L+1 + F + B + G ,
where F and B are found with the same expressions as the idle times for Case 1 of Forward phase (see Section B.1).

As no offloading is possible during backward propagation and no prefetching is possible during forward propagation, the idle time between the phases comes from the completion of both communication tasks, i.e. G = (∆ F L+2 + ∆ B L+2 )/β.

B.3 Backward phase

The situation in the case of backward is a bit more complex. We must indeed perform all the operations B j , . . . , B i+1 , with only I x i and δ j into memory. This is nevertheless enough since we can execute the whole forward chain F i+1 , . . . , F j from a i and then the whole chain B j , . . . , B i+1 using the values computed during the processing of the forward chain.

In general, due to memory constraints, it is not possible to perform F i+1 , . . . , F j B j , . . . , B i+1 in sequence. Since we assume that offloading cannot take place in the backward phase, we rely on rematerialization in order to save memory if needed. Our goal, given ∆ Fi and M Fi , is to find a valid schedule, that satisfies memory constraints and whose duration is minimal. One additional difficulty comes from prefetched data. Indeed, several valid rematerialization sequences for computing B j , . . . , B i+1 will differ both by their duration and by the amount of data that can be prefetched during the sequence. Indeed, for a given ∆ Bi , the sequences that enable to prefetch a lot of data during the computation of B j , . . . , B i+1 are preferable since they will induce a smaller value of ∆ Bj and therefore less memory pressure to perform B j , . . . B L .

Let C ∆ (i, j, A, ∆) denotes the optimal duration to execute the chain between layers i and j with A denoting the cumulative size of all activations already stored in GPU, including input I x i , and given minimal available memory M GPU -A -∆. In the beginning of the execution I x i and δ j are stored in the device memory. Moreover, I x i must stay in the memory until the end of the execution. Therefore, the first operation should be F i all or F i ck . Depending on the first operation, the optimal duration can be read in C F all ∆ (i, j, A, ∆) or C F ck ∆ (i, j, A, ∆) respectively. Case 1: F i all is the first operation Let us first notice that after F i all the activation āi must be kept in the memory until its associated backward operation (due to memory persistency). We must solve the subproblem of finding the optimal duration for the chain between layers i + 1 and j, which can be scheduled optimally with C ∆ (i + 1, j, A + āi , ∆ 1 ), according to our assumption. Once B i+1 is performed, then B i can be directly executed, as all necessary input are already in the device memory. Therefore, C F all ∆ (i, j, A, ∆) is given by C F all ∆ (i, j, A, ∆) = u Fi + C ∆ (i + 1, j, A + āi , ∆ 1 ) + u Bi + F i all + Bi .

The parameters for the recursive call of the dynamic programming are given by ∆ 1 = max{∆ -(u Bi + Bi )β, 0}.

We can also estimate the minimal available memory when F i all is executed, which is given by M GPU -A -∆ 2 where ∆ 2 = max{∆ 1 -C ∆ (i + 1, j, A + āi , ∆ 1 )β, 0}.

The idle times caused by prefetching can in turn be computed with the help of Lemma 1 as

Bi = max M i,j Bi -M GPU + A + ∆ β , 0 and F i all = max   M i,j F i all -M GPU + A + ∆ 2 β , 0   . M i,i F i all > M GPU -A or M i,i
Bi > M GPU -A means that there is not enough memory to perform this sequence, and as previously the dynamic programming cost should be set to ∞.

Case 2: F i ck is the first operation Let us suppose that F i ck is used to save I x i and consider the next value a i to be kept in memory. To compute this value, a sequence of F ∅ operations from layer i + 1 till layer i is performed. Due to memory persistency, after checkpointing a i we can find the optimal sequence from i + 1 to j by reading C ∆ (i + 1, j, A + a i , ∆ 2 ). Once B i +1 is performed, the chain from i to i can be scheduled using C ∆ (i, i , A, ∆ 1 ). Thus, C The parameters for the recursive call of the dynamic programming can be found using ∆ 1 = ∆ and ∆ 2 = max{∆ 1 -C ∆ (i, i , A, ∆ 1 )β, 0}.

Similarly to the previous case, we can estimate the minimal available memory when forwards from i to i are performed for the first time. This minimal memory is given by M GPU -A -∆ 3 , where

∆ 3 = max{∆ 2 -C ∆ (i + 1, j, A + a i , ∆ 2 )β, 0}.
The idle time caused by prefetching can be computed with the help of Lemma 1 as

F = max   max k=i,...,i (M i,j F k ∅ -β( i h=k+1 u F h )) -M GPU + A + ∆ 3 β , 0   .
If M i,j F k ∅ > M GPU -A for some k, i ≤ k < j, then this means that there is not enough memory to perform this sequence, and as previously the dynamic programming cost should be set to ∞.

Initialization: solution for a single layer Each recursive call to the dynamic programming is solving the problem for smaller sub-chains. The recursion stops when reaching the sub-chain consisting of only one layer. The schedule for one layer i is straightforward: perform F i all and B i , i.e. C ∆ (i, i, A, ∆) = u Fi + u Bi +

Figure 1 :F

 1 Figure 1: Data dependencies induced the training phase of Sequential Deep Neural Networks. Operation Input Output Time Memory overhead F all Forward and save all {a-1 } {a -1 , ā } u F o f {ā -1 } {ā -1 , ā }

Figure 3 :

 3 Figure 3: Simulation results for fixed memory ratio α = 0.2 and varying bandwidth.

  until it reaches m min , following the rule m(d) = m min + βd, where d denotes the distance to the end of prefetching. Therefore, for each operation o, d must satisfy m min + βd ≥ M S o and by construction d ≥ d o . Therefore, for each o, some idle time d -d o ≥ M S o β -d o -mmin β must take place before o. Since this constraint must be satisfied for all o, then = max max o∈S (M S o -βd o ) -m min β , 0 . B Proof of Theorem 1. Detailed dynamic program to compute the optimal sequence

F∆

  ck (i, j, A, ∆) is given byC F ck ∆ (i, j, A, ∆) = min i≤i <j i k=i u F k + C ∆ (i + 1, j, A + a i , ∆ 2 ) + C ∆ (i, i , A, ∆ 1 ) + F

  Simulation results for fixed bandwidth β = 12GB/s and varying memory ratio α.

		densenet 121	densenet 161	densenet 169	inception 200
	Overhead wrt sequential (for 12GB/s bandwidth)	1.0 1.1 1.2 1.3 1.0 1.1 1.2 1.3	resnet 101		resnet 152	resnet 200		resnet 50	Algorithm autocapper offload opportunist POFO rematerialization
		0.2 0.4 0.6	0.2 0.4 0.6	0.2 0.4 0.6	0.2 0.4 0.6
					Memory Ratio
	Overhead wrt sequential (for memory ratio α=0.2)	densenet 121 Figure 2: resnet 101 1.0 1.1 1.2 1.3 1.4 1.0 1.1 1.2 1.3 1.4	densenet 161 resnet 152	densenet 169 resnet 200	inception 200 resnet 50	Algorithm autocapper offload opportunist POFO rematerialization
		12 16	24	36 12 16	24	36 12 16	24	36 12 16	24	36
					Bandwidth (GB/s)

Table 2 :

 2 Values for the new state

https://pytorch.org/docs/stable/_modules/torch/utils/checkpoint.html

https://gitlab.inria.fr/hiepacs/rotor

https://www.deepspeed.ai

where, according to Lemma 1

Again, this solution will be infeasible if there is not enough memory to perform any operation, i.e. if

Bi > M GPU -A. In this case, the dynamic programming cost should be set ∞. Combining everything together: Therefore, C ∆ (i, j, A, ∆) can be computed as

B.4 Complexity

Finally, finding the optimal duration schedule for a chain of length L corresponds to computing OC x=0 (1, 0, 0, 0) using the dynamic programming presented above. The analysis of the complexity of the above dynamic programming can be decomposed in two parts. Let N GPU denote the number of discretized values for the memory M GPU . During the backward phase, the size of the state space is O(L 2 N 2 GPU ), and for case 2, computing a new value requires O(L) operations, which leads to O(L 3 N 2 GPU ) operations. In the forward phase, the size of the state space is O(LN 3 GPU ), and again computing a new value requires O(L) operations. This results in O(L 2 N 3 GPU ) operations for the forward phase. Therefore, the overall complexity is given by O(L 2 N 3 GPU + L 3 N 2 GPU ). In practice, we observe that for all the experiments presented in Section 5, discretizing the memory with N GPU = 50 values is enough since considering a finer discretization does not lead to any practical improvement.


Siavash Bakhtiarnia

Saeed Sheibani

Alain Billard

Hui Sun

Eric Aubry

Mohammad Arab

Pour Yazdi

Enhanced photocatalytic activity of sputter-deposited nanoporous BiVO4 thin films by controlling film thickness

Keywords: BiVO4, Thin film, Thickness, Sputtering, Photocatalyst

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Enhanced photocatalytic activity of sputter-deposited nanoporous BiVO4 thin films by controlling film thickness

Siavash Bakhtiarnia, Saeed Sheibani, Alain Billard, Hui Sun, Eric Aubry,

Mohammad Arab Pour Yazdi

To cite this version:

Siavash Bakhtiarnia, Saeed Sheibani, Alain Billard, Hui Sun, Eric Aubry, et al.. Enhanced photocatalytic activity of sputter-deposited nanoporous BiVO4 thin films by controlling film thickness. Journal of Alloys and Compounds, 2021, 879, pp.160463 (14). hal-03360033


Vinicius Ribeiro
email: vinicius.souza-ribeiro@loria.fr

Karyna Isaieva
email: karyna.isaieva@univ-lorraine.fr

Justine Leclere
email: justine.leclere@chu-reims.fr

Pierre-André Vuissoz
email: pa.vuissoz@chru-nancy.fr

Yves Laprie
email: yves.laprie@loria.fr

Towards the prediction of the vocal tract shape from the sequence of phonemes to be articulated

Keywords: phoneme-to-articulatory, speech production, neural networks

In this work, we address the prediction of speech articulators' temporal geometric position from the sequence of phonemes to be articulated. We start from a set of real-time MRI sequences uttered by a female French speaker. The contours of five articulators were tracked automatically in each of the frames in the MRI video. Then, we explore the capacity of a bidirectional GRU to correctly predict each articulator's shape and position given the sequence of phonemes and their duration. We propose a 5-fold cross-validation experiment to evaluate the generalization capacity of the model. In a second experiment, we evaluate our model's data efficiency by reducing training data. We evaluate the point-to-point Euclidean distance and the Pearson's correlations along time between the predicted and the target shapes. We also evaluate produced shapes of the critical articulators of specific phonemes. We show that our model can achieve good results with minimal data, producing very realistic vocal tract shapes.

Introduction

Continuous speech is a dynamic and non-stationary process that requires the interaction of several articulators. It is essentially the rapid transitions between vocal tract configurations that allows speech production [START_REF] Saha | Towards automatic speech identification from vocal tract shape dynamics in real-time mri[END_REF] and the articulation of phonemes is thus very context-dependent [START_REF] Öhman | Coarticulation in vcv utterances: Spectrographic measurements[END_REF]. A model of the vocal tract shape during speech is the subject of scientific interest for years. At first, Öhman [START_REF]Numerical model of coarticulation[END_REF] proposed a numerical model consistent with vocal tract measurements on X-ray images of Sweedish vowelconsonant-vowel utterances. Later, Maeda proposed an articulatory model [START_REF] Maeda | Un modèle articulatoire de la langue avec des composantes linéaires[END_REF] and then studied compensatory effects [START_REF]Compensatory articulation during speech: Evidence from the analysis and synthesis of vocal-tract shapes using an articulatory model[END_REF] showing how different articulatory pairs could produce the same phonetic pattern. Then, [START_REF] Potard | Incorporation of phonetic constraints in acoustic-to-articulatory inversion[END_REF] incorporated phonetic constraints derived from the common knowledge of inversion to the model.

Nowadays, deep learning dominates the field of speech processing, using extensive amounts of data, new learning schemes and computational power. Even though deep neural networks have an incredible impact in many areas, speech production still presents many challenges due to the diversity of physiological and/or physical phenomena involved and the high intra-and inter-speaker variability. [START_REF] Serrurier | Characterization of inter-speaker articulatory variability: a twolevel multi-speaker modelling approach based on mri data[END_REF] characterized speakerindependent articulatory strategies and inter-speaker variability for 11 French speakers uttering 62 vowels and consonants. The author's model is capable of explaining 66% to 69% of the variance. However, the research suggests that most of the variability is due to anatomical differences instead of speaking strategies.

Current research on the direct problem of phoneme-toarticulatory prediction, or its inverse form -acoustics-toarticulatory inversion -varies not only in the method itself but also in the source of articulatory data. The most common data acquisition modality is electromagnetic articulography (EMA), which has been extensively utilized in research [START_REF] Biasutto-Lervat | Phoneme-to-articulatory mapping using bidirectional gated rnn[END_REF][START_REF] Illa | Low resource acoustic-to-articulatory inversion using bi-directional long short term memory[END_REF][START_REF] Bozorg | Acoustic-to-articulatory inversion with deep autoregressive articulatory-wavenet[END_REF][START_REF] Shahrebabaki | Sequence-to-sequence articulatory inversion through time convolution of sub-band frequency signals[END_REF].

The disadvantage of EMA is the small number of sensor coils attached to the subject's articulators and the impossibility of measuring articulators located deeper in the vocal tract. Therefore medical imaging techniques are the most promising to account for less accessible articulators. MRI is a non-invasive imaging technique that provides a precise observation of all the vocal tract structures (except bones and teeth) without presenting any hazard for the patient. Current real-time MRI (RT-MRI) enable acquisitions with a high frame rate at a lower resolution [START_REF] Uecker | Real-time mri at a resolution of 20 ms[END_REF].

Csapó [START_REF] Csapó | Speaker dependent acoustic-to-articulatory inversion using real-time mri of the vocal tract[END_REF] explored RT-MRI for acoustic-to-articulatory inversion. Midsagittal RT-MRI images of the vocal tract were estimated using MGC-LSP spectral features as input and showed that the LSTMs are the most suitable for the task. Even though it is an excellent concept, the generated images contain several artifacts, and the produced shapes are not sufficiently realistic. We hypothesize that by targeting the complete MRI frame, the model misused the representative capacity to learn features that are not related to speech production.

In this work, we explore an alternative approach. We want to produce a direct mapping between the phonemes to be articulated and four vocal tract articulators -lower and upper lips, tongue, and soft palate -plus the pharyngeal wall which is useful to derive the nasopharyngeal port. Those five mentioned articulators are tracked in training images using a deep convolutional neural network approach (DCNN). The speech sound recorded with an optical microphone, and then denoised, enables ASR "forced" alignment and phonetic segmentation. Each frame is thus labelled with the corresponding articulated phoneme. Figure 1 presents a sample of the delineated contours automatically produced by our method [START_REF] Isaieva | Automatic tongue delineation from mri images with a convolutional neural network approach[END_REF]. This papers describes the approach used to predict the position of the articulators from the raw sequence of phonemes to be articulated.

Materials

Corpus

The dataset used corresponds to one adult female French native speaker. The RT-MRI sequences were acquired at the Centre Hospitalier Régional Universitaire de Nancy (CHRU Nancy), France. The recordings have a frame rate of 50 fps and image resolution of 136x136 pixels. The audio recordings are sampled at 16 000 Hz. We used Astali [START_REF]Astali[END_REF] to "forced" align the speech with the transcriptions and obtain the phonetic segmentation. An expert manually corrected the phonetic annotations.

The corpus contains 16 acquisitions, with a median acquisition time of 51.5 seconds, a minimum of 43 seconds and a maximum of 79 seconds. The acquisitions have a median of 7 sentences each, a minimum of 4 sentences and a maximum of 9 sentences. Some of the sentences are repeated during the acquisition. The sentences were selected to provide a phonetically balanced coverage of French. The corpus includes 107 spoken sentences, with 36 unique phonemes plus 4 tokens to represent silence, inaudible or non speech noises, hesitations, stop closures and burst releases.

Automatic tracking of articulators

To track the articulators, we trained a Mask R-CNN model [START_REF] He | Mask r-cnn[END_REF] to segment boundaries of the articulators. Then, we postprocess the probability map using the graph-based algorithm described by [START_REF] Isaieva | Automatic tongue delineation from mri images with a convolutional neural network approach[END_REF]. The algorithm searches for the two most distant points in the probability map and connects them using Dijkstra's shortest path algorithm. For the soft palate, our approach is subtly different. Because it is a thin structure, the algorithm does not perform correctly. We explored other approaches, e.g., active contours models [START_REF] Kass | Snakes: Active contour models[END_REF], but the results were not satisfactory. Instead, we modified the network's output such that the graph-based approach would generate the center line of the soft palate, which we use as target. In the end, the contours are regularized using b-splines to achieve a smoother curve. For the two cases in which the model could not extract the contour for the lower lip, we did it manually.

Generation of the vocal tract shape from

the sequence of phonemes

Architecture

We employed the architecture proposed by [START_REF] Biasutto | Modeling labial coarticulation with bidirectional gated recurrent networks and transfer learning[END_REF], with a few modifications. The network depicted in Figure 2 contains two parts. The first includes two layers of bidirectional gated recurrent units (GRU). Then, a linear layer with ReLU activation reshapes the GRU cell's output. The second part is the Articulator Predictor block, which comprises a sequence of layer normalization and linear layers with ReLU activation. The sigmoid activation in the end guarantees the outputed coordinates are between 0 and 1. Our architecture is composed by five Articulator Predictor blocks, one for each articulator. The input of the network is the phonetic segmentation, i.e. the phonemes to be articulated and their boundaries. In this current work there is no prediction of the sound duration from the sequence of phonemes to be articulated. The network produces the five contours for each 20 ms frame.

Evaluation

The models were trained by minimizing the Euclidean distance between the predicted spatial coordinates and the target positions. The loss function is given by

L(p, p) = 1 Nart × N samples N art i=1 N samples j=1 d(pij, pij) (1)
where p and p are the ground truth and the predicted points, Nart is the number of articulators, N samples is the number of predicted samples, and d is the Euclidean distance. In our case, Nart = 5 and N samples = 50. We also measure Pearson's correlation between the predicted and the target curves' trajectories along time.

It should be stressed the contours used as the ground truth are automatically delineated by the algorithm described in Section 2.2. Despite its very good robustness some contours might be erroneous.

Finally, we make a subjective analysis of the produced shapes. We evaluate specific phonemes that have critical articulators. Critical articulators are those resistant to context and that have a co-articulatory effect on neighboring phones [START_REF] Silva | Critical articulators identification from rt-mri of the vocal tract[END_REF]. For example, for the phonemes /p/ and /b/, it is mandatory to have the lips fully closed. The tongue tip must touch the teeth to produce the phonemes /l/ and /t/. For /k/, the tongue dorsum must approach the hard palate, and for nasal vowels, the soft palate must be open. We evaluate our model to check if these physical constraints are satisfied.

Experimental Setup

We propose two experiments. In the first, we separate three of the 16 acquisitions in our dataset for the test. From the remaining 13 acquisitions, we run a 5-fold cross-validation scheme, keeping 11 for training and two for validation. In the second, we train the model with 14, 11, 8, 6, 4, and 2 acquisitions, keeping validation and test sets fixed with one acquisition each.

We feed each sentence phoneme-by-phoneme to the model. To account for phoneme duration, we repeat the phonemes to match the number of frames of the RT-MRI films. Thus, each token in the input sentence lasts for 20 ms. The model predicts one shape per input token.

We trained our models for 300 epochs, with 20 epochs of patience for early stopping. We used the Adam optimizer [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF], with a learning rate of 1e-4, and we reduce the learning rate by a factor of 0.1 after ten epochs without improvements. We implemented the code using Pytorch [START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF].

Results

Table 1 presents the results in the 5-fold cross-validation framework. Figure 3 presents our results in the data efficiency experiment. Figure 5 depicted figures present the predicted contours after b-spline regularization, which we apply only for visualization.

In our second experiment, i.e., data-efficiency evaluation, the larger experiment -with more acquisitions during training -fully contains the data from all of the smaller ones. Thus, the number of training sentences increases with the number of acquisitions and phonemes present in the smaller experiments are present in the larger ones as well.

Table 1: Mean Euclidean distance, X, and Y correlations for each articulator in the 5-fold cross-validation scheme. We omit the mean Euclidean distances' standard deviations to save space, but they are approximately 0.001.

Articulator MED ρx ρy

Tongue 0.016 0.689 ± 0.065 0.637 ± 0.050 Pharynx 0.007 0.327 ± 0.099 0.396 ± 0.047 Soft palate 0.012 0.407 ± 0.031 0.403 ± 0.010 Lower lip 0.011 0.668 ± 0.050 0.714 ± 0.028 Upper lip 0.008 0.535 ± 0.052 0.512 ± 0.034

Discussion

Table 1 and Figure 5 show that the model can produce highquality shapes, with tiny errors and good correlations between the predicted and target trajectories. Even though the lowest errors occur for the upper lip, the best x-and y-correlations occur for the tongue and the lower lip. The tongue is the largest considered articulator, with the highest number of degrees of freedom. Thus it is expected that it accounts for the most significant errors. The pharynx accounts for the lowest correlations, which is expected since it is a fixed articulator, with almost no degrees of freedom. Figure 4 shows that the model fails to reconstruct the sublingual cavity when it is present in the ground truth contours. We hypothesize that it is an error propagated from the tracking algorithm which sometimes fails to identify the sublingual cavity. Usually, this is a fuzzy region in the low-resolution image and the probability map is not sufficiently clear.

On the data efficiency experiments, we evaluate the model's capacity of learning from fewer data. Our results show that our model is highly data-efficient. Even though we reach the best results using all the acquisitions during training, the marginal gain with more than six acquisitions is less expressive. This behavior is observable for both Euclidean distance and correlations.

From the critical articulators' perspective, our model can learn some specific physical constraints. As discussed before, lips should be closed for the production of /p/ and /b/. Figure 5 shows that this condition is satisfied for /p/ but not for /b/. Actually there is a closing movement but it is not complete. One reason is probably the fact that the closure for /b/ is shorter than /p/, thus corresponds to very few frames with a complete closure, with the consequence of not training the model with enough relevant data.

The evaluation is trickier for the phonemes /k/, /t/, and nasal vowels because they depend on the interaction with fixed structures. For /t/, we see that the model projects the tongue tip forward, suggesting contact with (virtual) teeth. However, for some cases in the dataset, we expected a more prominent tongue tip. For /k/, a similar effect occurs with the tongue dorsum, which the model projects up, closing the gap with a (virtual) hard palate.

When we compare the last row of Figure 5 with the previous ones, which refers to the /ã/, we observe that the soft palate opens for nasal vowels, and closes for non-nasal phonemes.

When we compare the predicted curves with the ground truths, we see that the critical articulator constraints are mostly satisfied in our model's predictions.

Conclusions

In this paper, we present a novel approach for predicting the vocal tract shape from the sequence of phonemes to be articulated. To the best of our knowledge, this work is the first to produce the main articulator contours for continuous speech. Our model can yield promising results for a single speaker with minimal data even if one expects better results with a bigger corpus. In the subjective evaluation, critical articulators' physical constraints are satisfied, which shows that our results are realistic.

However, there are still limitations. First, this research focuses on articulators in the upper part of the vocal tract. An obvious continuation is the addition of other structures, i.e. the epiglottis and larynx, which are essential to get the complete geometry so as to synthesize a speech signal by means of the numerical simulations of the aero-acoustic equations.

Second, we rely on automatic contour tracking which produces reference contours used as ground truth. The presented model has a performance upper bound linked on the capacity of producing an accurate and relevant ground truth.

Last, we work with a single speaker. Our input does not encode any speaker-specific feature, and it only uses the sequence of phonemes to be articulated and their expected duration. Thus, learning from several speakers is a challenge and a sufficient amount of data is a requirement to marginalize anatomical variance.

Figure 1 :

 1 Figure 1: Sample of MRI image with lower lip (green), upper lip (magenta), tongue (orange), pharynx (dark yellow), and soft palate (blue) and wireframe representation contours. Contours were automatically extracted from MR images.

Figure 2 :

 2 Figure 2:The proposed network architecture. We present the recurrent encoding path on the and the articulator predictor decoder on the right. Each RNN cell output is reshaped and fed into the articulator predictor. We draw only the first cell for simplicity.

Figure 3 :Figure 4 :

 34 Figure 3: Results from the data efficiency experiment. On the left, the mean Euclidean distance for each articulator. On the right, the X (dashed line) and Y (dotted line) correlations for each articulator. /u/ /@/

  presents samples of critical articulators. The

	5x	Articulator Predictor
	BiGRU BiGRU	Linear ReLU Latent Space		... ...	Latent Space	LayerNorm	Linear	ReLU	LayerNorm	Linear	ReLU	LayerNorm	Linear Linear	X Coords Y Coords
		#	t	u	#

  Figure 5: Samples of our model predictions for the phonemes /b/, /t/, /p/, /k/, and /ã/, in order of appearance from top to bottom. The solid lines represent model predictions. The red dotted line represents the ground truth. Each row represents one phoneme occurrence, and the columns represent the time steps (20 ms).

	/b/	/b/	/b/
	/t/	/t/	/t/	/t/
	/k/	/k/	/k/	/k/
	/p/	/p/	/p/	/p/	/p/
	/a~/	/a~/	/a~/	/a~/	/a~/


EDITORIAL

Pioneers of Plankton Research: Marie Lebour (1876Lebour ( -1971) ) Marie Lebour made major contributions to plankton research. Her books on dinoflagellates and diatoms were the first in the English language and remain today valuable reference works. Her studies of larval fish of the plankton, both as predators in the plankton and the prey of other planktonic organisms, have been cited hundreds of times (e.g. [START_REF] Lebour | The food of plankton organisms[END_REF], "The Food of Plankton Organisms"). However, Marie Lebour is today perhaps best known from her pioneering work on the identity and developmental sequences of many of the striking larval forms of crustaceans and mollusks we find in the plankton. The following brief resumé of Marie Lebour's long and remarkably productive career is based largely on the account of her life given in Russell (1972), which includes a complete bibliography of Lebour, and the account given in [START_REF] Rice | Marie Victoire Lebour (1876-1971)[END_REF].

Marie Lebour appears to have come to plankton studies via a very unusual route. As a young naturalist, her first published study was on land and freshwater mollusks [START_REF] Lebour | Land and freshwater Mollusca in Northumberland, collected chiefly at Corbridge-on-Tyne and on the lower Tweed[END_REF]). She published 3 other studies on marine mollusks before finishing her undergraduate studies in 1904 at the University of Durham. By 1906 she was employed as a 'junior demonstrator' in the Department of Zoology at the University of Leeds. She remained at the University of Leeds for several years, publishing on near-shore mollusks and their parasitic trematodes. While she concentrated mostly on trematodes of mollusks, she also worked on the trematodes of ostrocod crustaceans, fish, and even a snake. In this period, her work on trematodes, and their development in their hosts, appear to be her first studies in developmental biology. Decades later, she would become known as an expert in the development of a very wide variety of marine planktonic larvae.

Her studies of plankton began only in 1915, at age 39, when she went to Plymouth Marine Laboratory (PML). Oddly enough, the move to the PML was not one she sought. Marie Lebour was 'loaned' by the University of Leeds to PML supposedly for only a year, where at the time of WW1, a large part of the staff had been called to military service. However, PML apparently suited her quite well. In 1917 she was offered a permanent position at the University of Leeds but chose to stay at the PML. She would remain at the PML for the rest of her career, working in the laboratory until 1964, well past her official retirement in 1946. In fact, post-retirement she published 17 papers. Her last, and 156th paper appeared in 1959 [START_REF] Lebour | The larval decapod Crustacea of tropical West Africa[END_REF], when she was 83 years old! One can roughly parse the career of Marie Lebour as a plankton specialist at PML into three more or less distinct periods, which is to say with some temporal overlaps. In the first period were her studies of microplankton, chiefly dinoflagellates and diatoms of the phytoplankton. In the second period, her major works were investigations into the biology of larval and young post-larval fish, both as predators, (studying what they ate from fish gut contents), and as prey of other planktonic organisms in laboratory studies with invertebrate predators. During the third and last period the studies of Marie Lebour largely concerned documentation of the development stages of a wide variety of the meroplanktic larvae of both molluscan and crustacean taxa. Marie Lebour's first published work on marine plankton was her now classic study documenting the developmental stages of the copepod Calanus finmarchicus [START_REF] Lebour | Stages in the life history of Calanus finmarchicus (Gunnerus), experimentally reared by Mr L. R. Crawshay in the Plymouth Laboratory[END_REF]). The copepod had been brought into culture by L.R. Crawshay. He had painfully worked out culture methods over several years at PML. When Crawshay was called to military service, he gave his cultures to E.J. Allen, the director of PML at the time, and they supplied the material used by Lebour for her first publication from PML, signed as "Assistant Lecturer in Zoology, Leeds University. Temporary Naturalist in the Plymouth Laboratory." in 1916 she also began her work based on organisms collected in plankton net tows, working out developmental stages of a crab parasitic in medusa. However, she did not abandon her work on trematodes as she also described trematodes she found as larvae in medusa.

In just the first few years at PML, Marie Lebour attained a considerable expertise in the study of plankton. By 1917 Marie Lebour had become the British expert on the dinoflagellates of the marine plankton, describing several new species in her first work on dinoflagellates [START_REF] Lebour | The peridiniales of Plymouth Sound from the region beyond the Breakwater[END_REF]. A few years later she published "The Dinoflagellates of the Northern Seas" (Lebour 1925a), her first book on phytoplankton and the first book in English on dinoflagellates [START_REF] Rice | Marie Victoire Lebour (1876-1971)[END_REF]. It was illustrated with her detailed drawings of cells from life (e.g., Fig. 2) and was very favorably reviewed both in the United Kingdom (Anon 1926a) and in the United States of America (Anon 1926b). Marie Lebour's contributions to our knowledge of dinoflagellates is evidenced by the fact that several genera and species of dinoflagellates have been named for her. Among those currently recognized as valid are the genera Lebouraia and Lebouridinium and the species Cochlodinium lebourae and Polykrikos lebourae. The fact that phagotrophic species of dinoflagellates (Cochlodinium and Polykrikos) are named for her is especially apt. She was fond of insisting that dinoflagellates are not typical plant-like phytoplankton but rather at borderland separating plants and animals with some forms feeding and other relying on photosynthesis (e.g., Lebour 1928). Marie Lebour was an early advocate of the fact that within a taxon of marine plankton, a variety of trophic modes can be found. It was also during those early years at PML, following her work on dinoflagellates, when Marie Lebour began her studies of larval and young post-larval fish as both predators and prey of planktonic organisms. The first reports in 1918, signed as "Naturalist at the Plymouth Laboratory" (no longer 'Temporary Naturalist'), were on the prey found in the alimentary canal of the in a large variety of species. She found that certain fish appeared to prefer certain prey (e.g., [START_REF] Lebour | The food of post-larval fish[END_REF]. She continued her studies on the food of larval fish, based on field-caught material, over the next few years finishing with studies of the larvae of commercially important species such as Herring [START_REF] Lebour | The larval and pest-larval stages of the pilchard, sprat and herring from Plymouth district[END_REF]. Thereafter, she turned her attention to feeding of other organisms of the plankton, primarily coelenterates, on fish larvae. These studies relied on laboratory observations and were made on field-caught individuals maintained in the laboratory in vessels she developed, Plunger Jars, adapted from designs originally developed at PML to maintain medusae. The Plunger Jars kept the organisms in suspension. She documented predation on fish larvae by a surprisingly wide variety of invertebrates ranging from medusa, ctenophores (fig. 3) to arrow worms and annelids. Her last publication on holoplanktic organisms was her 1930 book on diatoms, "The Planktonic Diatoms of Northern Seas" [START_REF] Lebour | The Planktonic Diatoms of Northern Seas. Ray Society[END_REF]. Subsequently her studies focused on the meroplankton, specifically the development of the planktonic larvae of benthic taxa. From the 1930's on, Marie Lebour 's work largely concerned documentation of the developmental stages of crustacean larvae. She produced a large number of studies such that crustacean biologists often consider her to be one of their own, ignoring her previous and other work [START_REF] Rice | Marie Victoire Lebour (1876-1971)[END_REF]). However, during this third part of her career, she also worked on the planktonic forms and stages of mollusks, the group of organisms with which she began from her days as a young naturalist. Her renown among crustacean biologists is due to the fact that she was one of the first to successfully rear crustacean larvae and her studies, unlike those of her predecessors which were largely morphological, examined questions of nutrition, ecology, as well as patterns of temporal and spatial distributions [START_REF] Anger | The Biology of Decapod Crustacean Larvae[END_REF]. Many of her studies included illustrations in color, showing the natural colors of the larval forms, such as those shown in Figure 4. She apparently studied art as an undergraduate and her artistic talents were well recognized (Russell 1972) In concluding our review of the work of Marie Lebour, it would be remiss to omit mention of her work in popularizing the world of plankton. She contributed several articles over the years to the popular science magazine Science Progress. One of her articles told of the story of her exploiting the fact that her brother-inlaw was traveling to the then Ceylon (now Sri Lanka) on a cruise ship equipped with a running seawater bath. The article was titled "Through the Bath Tap (an account of samples of minute planktonic organisms collected on board a P & O liner during a voyage from England to Ceylon)". She supplied her brother-in-law with a bag of silk fashioned to fit over the bath faucet and tubes with a small amount of formalin in each. Periodically he tied the bag on the faucet and let it run for about an hour. He then noted the position of the ship, and emptied the netting into a tube. From her article, the "bath tap plankton net" and some of the catch are shown in Figure 5. Thus, the readers of Science Progress were shown the variety of life-forms found in the plankton of tropical seas. Although not identified in her article, the illustrations were detailed enough to allow us today to easily identify many of the species shown. She also wrote several articles for a magazine for science teachers, School Science Review, for example "Plankton or the Floating Life of the Sea" [START_REF] Lebour | Plankton, or the floating life in the sea[END_REF]). Marie Lebour was not only a pioneer of plankton research but also a public proponent of the importance and beauty of the plankton.

Fig. 1 .

 1 Fig. 1. Marie Lebour ca. 1916, at about age 40. Image from Russell's 1972 obituary notice published in the Journal of the Marine Biological Association.

Fig. 2 .

 2 Fig. 2. Plate 12 from Lebour's "The Dinoflagellates of the Northern Seas" (Lebour 1925a).

Fig. 3 .

 3 Fig. 3. From Lebour's 1925b study "Young anglers in captivity and some of their enemies, a study in a plunger jar" in the Journal of the Marine Biological Association (UK), reproduced by permission.

Fig 4 .

 4 Fig 4. Color plates of Marie Lebour from her 1931 article on larval Brachyura (left panel) and her 1935 article on Echinospira larvae (right panel). Both article appeared in the Proceedings of the Zoological Society (London), reproduced by permission.

Fig. 5 .

 5 Fig. 5. Illustrations from Marie Lebour's 1932 article "Through the Bath Tap (an account of samples of minute planktonic organisms collected on board a P & O liner during a voyage from England to Ceylon)" which appeared in the magazine Science Progress.