1. Introduction
In 2016, the US Department of Energy (DOE) established the Exascale Computing Project (ECP, https://www.exascaleproject.org) to accelerate delivery of a capable exascale computing system. In the ECP project, significant effort is being invested to develop highly scalable numerical libraries and highfidelity modeling and simulation codes across the spectrum of science and engineering domains and disciplines. Given a myriad of diverse computer architectures, achieving optimal performance and performance portability of all these codes has become an increasing challenge. Each of these codes has a number of tuning parameters that may be difficult to choose to optimize performance (or any quantitative metric). We are working on building an adaptive execution framework applicable across ECP applications. Adaptive execution strategies are concerned with picking the right set of parameters to optimally solve a particular problem on a given architecture. Our interrelated goals are to make it easy to do the following: 1) Support a variety of tuning metrics, including time, memory usage, accuracy, or hybrids like minimizing time given a memory constraint; 2) Compose large applications out of multiple code bases, including the ability to optimize the overall application, perhaps by changing data structures or the number of processors in between calls to different codes, or using different hardware resources; 3) Support archiving and reusing tuning data from multiple executions to allow tuning to improve over time; 4) Allow incorporation of new optimization techniques to improve the tuning process.
Since most of the exascale applications involve expensive “function evaluations”, requiring either long runtime or many hardware resources (e.g., core count), the bruteforce “grid search” approach to finding optimal parameters is infeasible. Our approach is to build an automatic optimization routine to choose these parameters, using the application code as a “black box”, and perform a small number of carefully chosen runs with different values of tuning parameters in order to find the best ones. These parameters could be internal to one routine, or involve multiple routines that must be “cotuned” because they share resources. If a routine is independent of, or has a low impact on, the rest of the application in terms of the cost metric to be optimized, this routine can be tuned separately, providing the user can execute it in a stand alone blackbox fashion. Otherwise, the full exascale application needs to be run in order for the autotuning to proceed.
As subpart of the ECP project, we have built a common interface for autotuning with the following features: 1) Input data and parameters, or any statistics describing them, are exposed by the application to the optimizer. These could include a pointer to prior information from a database, which we may also update with new data to improve the future adaptation layers; 2) Any resource constraints (e.g., core count or memory) are also exposed to the optimizer. These could come from the application or a higherlevel runtime system; 3) Optimization metrics and constraints (e.g., minimize runtime subject to a memory constraint) are also stated. Note that this data could be the union of data provided by many components of a larger application.
The core idea, and the main contribution, of this paper is to broaden the scope of autotuning, contrasting with that of stateoftheart approaches. Specifically, instead of tuning a given application on a particular target problem instance, as has been done so far, we aim to tune that application for a myriad of problems it aims to solve. The proposed realization of this idea is achieved through the use of multitask and transfer learning. Experimental results show a average improvement of the application runtime over stateoftheart autotuning approaches, together with competitive (runfree) autotuning model predictions. Moreover, we show that when the budget (number of allowed runs of the application) for autotuning is low, due to a high cost of the target application (e.g. exascale applications), our methods are more suitable than stateoftheart methods that are not necessarily designed for this use case.
The remainder of this paper is organized as follows. Section 2 describes the related work in autotuning and blackbox optimization for autotuning. Section 3 presents the multioutput autotuning framework together with the notations used in this paper. Section 4 and 5 explain the new methodology that relies on multitask learning and transfer learning, respectively, for the purpose of autotuning. Section 6 shows comparative results with the OpenTuner and HpBandSter stateoftheart autotuners, together with the new results which are not possible to get through standard approaches alone. Finally, section 7 summarises this work and discusses its perspectives and future directions.
2. Background
For a thorough survey of autotuning, multitask learning and transfer learning, we refer the reader to [8423171], [zy18] and [py10], respectively. Section 2.1 describes blackbox optimization as a mathematical framework for autotuning. Blackbox optimization techniques can be grouped in two categories. Modelfree or nonBayesian optimization methods try to find an optimum through a deterministic and / or stochastic decision process. Section 2.2 gives a brief overview of such methods, which are used in OpenTuner [opentuner] and HbBandSter [bohb], two stateoftheart autotuners that we compare our results with. Modelbased or Bayesian optimization methods build and optimize a surrogate model of the real objective function in order to find its optimum. Section 2.3 explains the principle of such methods, backbone of the machine learning methods of significance to our work, which are described later on in Section 3, along with our proposed methodology. Optimization methods that take into account applicationspecific knowledge are left aside to keep the focus on general purpose approaches instead.
2.1. Blackbox optimization for autotuning
The fundamental aspect of autotuning is optimization, i.e., finding a configuration of the parameters of an application that makes it solve a given input problem optimally. The nature of autotuning makes this optimization problem lie within the family of blackbox optimization problems, which are among the hardest to solve. Indeed, an expensive run of the application is necessary in order to get the value of the objective function to optimize (e.g. runtime, energy) for a given combination of parameters. Moreover, even in the presence of approximate (coarse) analytical models of the application (e.g. flop count, memory consumption), no precise enough formulation generally exists to predict more detailed phenomena (e.g. memory hierarchy communications, network contention) influencing the application behaviour. Furthermore, the human cost of deriving such models (if even possible to derive) is too expensive to be practical. Similarly, no analytical information on the gradient of the objective function is available, but only numerical approximations are, through sample evaluations of the objective function.
2.2. NonBayesian optimization
The simplest blackbox optimization methods, that are usually tried first before resorting to more advanced methods, are known as: The deterministic exhaustive search (and its variant grid search
), which tries all (or subset of all, respectively) possible combinations of all possible values of the parameters and selects the best performing one. Their drawback is that they quickly become intractable when the number of parameters increases, due to the curse of dimensionality
[Bellman:1957]; The stochastic random search, which selects randomly the value of each parameter in order to generate candidate solutions, then selects the best performing candidate.Two main families of modelfree optimization approaches exist. The global approaches explore the whole search space and try to find a balance between exploration of new regions of the search space and exploitation of the data gathered to give more attention to the promising regions of space. Examples of such methods are Simulated Annealing (SA) [sa], Genetic Algorithms (GA) [ga] and Particle Swarm Optimization (PSO) [pso]. In contrast, the local approaches try to improve upon previous solutions by exploring their neighboring region only until converging to a local minimum. Examples of such approaches are Nelder–Mead simplex [NeldMead65] and Orthogonal Search [clp11].
OpenTuner [opentuner] and HbBandSter [bohb]
are two stateoftheart general purpose autotuning frameworks. They rely on metaheuristics to solve a multiarmed bandit problem
[kv87] where application runtime is the resource to be allocated (in our case). OpenTuner allocates and distributes application runtime over a variety of optimization methods (mentioned in the previous paragraphs) in such a way as to adaptively select the best performing method, which is selected to solve the autotuning optimization problem. HpBandSter’s algorithm mixes a Bayesian Optimization method (see following section) with that of Hyperband [hyperband]. Hyperband, as an earlystopping method, allocates runtime uniformly over randomly sampled configurations, keeping only the half best performing configurations at each iteration, and extending the corresponding runtime per remaining configuration, until a single (best) configuration remains.Several other nonBayesian blackbox optimization packages for autotuning exist in the literature. Particularly, SuRf [surf]
uses random forests to model the performance of an application and find its optimum. One of its main strengths is its ability to handle categorical parameters (choices) in an elegant way.
2.3. Bayesian optimization
Bayesian optimization [sswaf16], also known as response surface methodology, relies on a surrogate model of the real objective function to optimize. A prior belief model representing the assumptions on the objective function is chosen, and a posterior is built from it so as to maximize the likelihood of some sample data of the objective function to be a realization of that model. Instead of directly optimizing the true objective function, the model is optimized instead, as it is much cheaper to evaluate, and iteratively updated until convergence to an optimum. The Efficient Global Optimization (EGO) algorithm [jsw98] is a classical Bayesian optimization algorithm. In order to balance between exploration and exploitation (global and local search behaviors), EGO tries to optimize a certain acquisition function (e.g. Expected Improvement (EI)) instead of the model itself. This metric considers both the value of the model at a given point in the search space together with the confidence in the model prediction at that location. For instance, if the model predicts a good value of the objective function with a high confidence at a given location, while a worse value at another location but with a lower confidence, it might be worth exploring this second location as an optimum might be located nearby. After finding a location in the search space that optimizes the acquisition function, this location is evaluated through the expensive blackbox objective function and the corresponding value is used to update the surrogate model. EGO iterates the process until the EI reaches a certain threshold, at which point the algorithm is considered to have converged to an optimum.
3. MultiOutput AutoTuning
Throughout the paper, we refer to input problem as an input of the target application to be tuned, that is. Also, we refer to task as the problem of tuning the parameters of the target application given a specific input problem.
Starting from the observation that tuning an application for a specific input problem makes that application efficient on solving that particular problem but without any guarantee on its behaviour on other unencountered problems, our aim in this paper is to tune an application for any problem it might encounter. This goal being potentially unrealistic in general, given the infinity of input problems that might exist, our aim is to find, in a reasonable amount of time, good enough combinations of parameters of the application for any given input problem.
The guiding ideas behind our work are twofold.

Firstly, we tune the application on a finite set of well chosen input problems. We believe that tuning the application on each problem independently from the others is less efficient than tuning all of them simultaneously, benefiting from the gathered knowledge on all input problems to speed up the whole tuning process This is critical especially in an exascale setting, where every run of the application is extremely expensive. We rely on the multitask learning framework to this end as detailed in Section 4.

Secondly, we derive from the optimal parameters found on the initial set of problems, together with the model built on these problems, approximations to the optimal parameters of any other (nontuned) input problem. We rely on the transfer learning framework to do so as explained in Section 5.
Let us define autotuning in the multioutput setting. Throughout this paper, the QR factorization routine of ScaLAPACK [slug], denoted as , is used as a guiding example of application to be tuned. It takes as input a matrix and computes its QR factorization, producing an upper triangular matrix and a matrix . Moreover, the notations used in this paper are summarized in Table 1.
Symbol  Interpretation  
General notations  
task space  
input space (parameter configurations)  
output space (e.g., runtime)  
dimension of  
dimension of  
dimension of  
number of tasks  
number of samples per task  
matrix of tasks selected from sampling  
matrix of samples (parameters)  
vector of results (runtime)  
ScaLAPACK related notations  
Task 
number of matrix rows  
number of matrix columns  
number of compute nodes  
number of cores per compute node  
Input 
row block size  
column block size  
number of MPI processes  
number of threads (used in BLAS)  
number of row processes  
number of column processes 
Let us define , the Task Space, as the space of all the input problems that the application may encounter. Let us make the simple (non restrictive) assumption that may be characterised as a finite dimensional space of dimension . This means that it is possible to identify any input problem with a finite number of features. Relating to , is a space of dimension . The corresponding features are: and , respectively the number of rows and columns of a matrix to be factorized; and , respectively the number of compute nodes together with the number of processor cores per compute node, characterizing the architecture of the machine on which the application runs. Several application problems are amenable to such a formalism, potentially after some approximations are made. If this finite dimension task space assumption does not hold, the multitask learning methods in Section 4 remain valid, but some of the transfer learning methods in Section 5 would require additional attention that we leave as a future work. An example of application where this assumption does not hold is the case of SuperLU [lidemmel03, Li05], where the task space represents the space of all possible sparse matrices, which cannot be characterized by a finite set of parameters. The underlying assumption that we make is that the objective function to optimize is somehow continuous in the task space, as we expect it to be similar for similar tasks (for the same parameters values).
Let us define , the Input Space, or space of the parameters to be optimized, as , with the number of parameters. In practice, inputs can be either real, integer or categorical (e.g. a list of algorithms (which can map to an interval )). In our work, the last two cases are translated internally to the real case. Every point in can be referred to as a parameter configuration. Moreover, most applications have constraints on their parameters as not all possible combinations of parameters (point in the input space) are valid. Relating to , is a space of dimension . The corresponding parameters are: and , the row and column block sizes, respectively; and , the number of MPI processes and BLAS threads, respectively, and and , the number of row processes and column processes, respectively. Moreover, two constraints on these parameters exist: () ; () . These constraints reduce the effective dimension of the input space to , as and can be deduced from the other parameters.
Let us define , the Output Space, as the space of the results of the evaluation of the optimization objective function, for a given task and on a given parameter configuration. It is a single dimensional space (e.g., computation time, memory consumption, energy …), the case of multidimensional spaces (corresponding to multiobjective optimization) being left for future work.
We denote by and the values of the objective function measure and model prediction , respectively, for a task and for a parameter configurations . As is the case in Bayesian optimization, the model is optimized instead of the measured value , while the optimum found is hoped to be that of . In this setting, given a task , the autotuning goal is to find:
(1) 
The autotuning stopping criteria is classically defined as one or a combination of the following: a maximum number of evaluations of the objective function (runs of the application); a maximum wallclock time; threshold on a relevant measure of the quality of the solution provided by the application (e.g., numerical accuracy, energy consumption). In order to be able to fairly compare different autotuning algorithms, we chose to fix a maximum number of runs of the application as the stopping criteria. The total runtime spent in the application is also measured and reported.
4. Multitask learning
In the field of machine learning, multitask learning consists of learning several tasks simultaneously while sharing common knowledge between them in order to improve the prediction accuracy of each task and / or speed up the training process. In this section, we propose the Multitask Learning Autotuning method, . The methodology followed in resembles that of singletask Bayesian optimization methods, as it adapts the EGO algorithm to the multioutput setting. It is composed of three main phases: sampling phase (Section 4.1); modeling phase (Section 4.2); and optimization phase (Section 4.3). The second and third phases are repeated until convergence or stopping of the tuning process.
4.1. Sampling phase
While a single sampling step is needed in a singletask Bayesian optimization scheme, two sampling steps are needed in .
The objective of the first sampling step is to select a set of tasks . The goal when selecting this set is to get a representative sample of the variety of problems that the application may encounter, rather than focusing on a specific type of problem. Given the freedom in the selection of the tasks together with the existence of a space of tasks , we chose a space filling sampling in to select the tasks. Such samplings are widely used in the field of Design Of Experiments (DOE) [doe]. Particularly, we chose a Latin Hypercube Sampling (LHS) [lhs] in our method. Such samplings try to cover the whole search space uniformly, such that no region of the space is oversampled while another one is left undersampled. Several offtheshelf software packages exist that implement several different types of samplings (including LHS).
The objective of the second sampling step is to select an initial sampling of parameter configurations for every task . For task , its initial sampling consists of parameter configurations . Two cases arise in the multitask framework: the isotropic and heterotropic cases. The isotropic case arises when all the tasks share the same sampling in the input space, while in the heterotropic case, different tasks do not necessarily share the same samples. In the isotropic case, the advantage of a multioutput regression is the sharing of information for the optimization of the hyperparameters of the model governing the tasks. In the heterotropic case, however, more knowledge can be shared as insights on the true cost of a task on an unknown configuration can be learned from a similar task with a sample at that location. Moreover, in reallife applications, given the existence of constraints on the parameters, not all parameter configurations are feasible for all tasks simultaneously; a configuration may be valid for a subset of tasks, but violate the constraints on another subset. Thus, we chose to generate the initial sampling in a heterotropic way by generating the as independent LHS.
Every sample is evaluated through a run of the application, whose result is represented as . The set represent the results of all these evaluations, , where every represents the results corresponding to task , .
Given the application constraints, a generic sampling technique might fail, both for the selection of the task samples and for the selection of their corresponding input samples. Such a case arise frequently when the total number of combinations of parameter values is of the order of thousands of billions (simply because of the curse of dimensionality) while the number of valid parameter configurations that respect the constraints is only of the order of hundreds of thousands. An example is the tuning of matrix multiplication on GPU in MAGMA [ahkld15, magma].
In such a case, either specific knowledge of the application should be used to design a space filling sampling, or, a simple remedy is to generate several unconstrained samplings and randomly select candidates from the union of sets of valid configurations drawn from every sampling.
4.2. Modeling phase
Once and are selected and evaluated, the core of is its modeling phase, which consists in training a model of the blackbox objective function relative to tasks . However, instead of building a separate model for every task, as is customarily the case in a regular singletask Bayesian optimization scheme, the challenge in is to derive a single model that incorporates them all, sharing the knowledge between them to be able to better predict them all. While Gaussian Processes (GPs) are often used in the modeling for single task tuning (e.g. [spearmint]), we propose to rely in on the generalization of GPs to the multioutput setting, known as the MultiOutput Gaussian Process framework [bkw08]. Particularly, we chose to use the Linear Coregionalization Model (LCM) [nla.catvn1116468, CIS7228], the choice of which is driven by its generality, flexibility and modeling power, albeit its modeling cost, compared to the plethora of models in the literature that are derived from it as special and constrained cases [bkw08].
A similar but orthogonal line of work that also relies on an agglomerate model of tasks is presented in [fsk07]. Indeed, in the context of multifidelity optimization, where one seeks the optimum of an expensive task while a hierarchy of highly correlated and cheaper tasks exist, the authors show that the use of coKriging (which can be formulated as a special case of LCM) leads to a more accurate model, and subsequently, to a faster convergence to a global optimum of the expensive task. However, instead of a (vertical) hierarchy of correlated tasks of increasing cost, our problem corresponds to a (horizontal) relationship between correlated tasks of the same (or similar) cost.
4.2.1. Gaussian Processes
We provide here a brief explanation of Gaussian processes. We invite the reader to consult [rw05]
for a detailed description. A Gaussian process is the generalization of a multivariate normal distribution to an infinite number of random variables. It is a stochastic process where every finite subset of variables follows a multivariate normal distribution. While other regression methods set a prior on the function to be predicted and try to learn the parameters of such a function, GPs set a prior on some characteristics of the functions (e.g. smoothness) and try to learn the functions themselves, allowing for the expressiveness of a much richer variety of functions. A GP is completely specified by its mean function
and by its covariance function . A function following such a GP is written as:(2) 
where:
(3) 
(4) 
In most practical scenarios, is taken to be the null function and all the modeling is done through , known as the kernel function. While a variety of kernel functions exist in the literature, the choice of which is dependant on the data to be modeled, a very popular generic choice is the exponential quadratic kernel given by:
(5) 
where
(variance) and
(length scales) are hyperparameters of the kernel governing its behavior. These are learned by optimizing the loglikelihood of the samples with values on the GP, the loglikelihood being described as:(6) 
where is a regularization term, and is the covariance matrix whose elements are generated from the kernel .
4.2.2. Linear Coregionalization Model
The key to LCM is the construction of an approximation of the covariance between the different outputs of the model (model of every ).
In this method, the relations between outputs are expressed as linear combinations of independent latent random functions
(7) 
where () are hyperparameters to be learned, and are the latent functions, whose hyperparameters need to be learned, as well.
The independence hypothesis is important as it allows us to compute the covariance between the outputs only through the autocovariance of the latent functions themselves, as it implies:
(8) 
The covariance of every latent function is assumed to be generated from a kernel function:
(9) 
This kernel can take and to be vectors as input, in which case, is a scalar. However, this kernel can also consider them to be vectors of vectors (i.e.: matrices). In such a case, is a matrix, instead of a scalar, where the entry corresponds to the evaluation of the kernel on the vector of and vector of .
The covariance structure between two outputs can now be expressed as:
(10) 
which, thanks to Equation (8) simplifies to:
(11) 
The covariance matrix between all the tasks on inputs and can now be expressed through the kernel
(12) 
where is a Kronecker product, is the covariance function of the latent function,
is the identity matrix of size
, is a diagonal matrix, whose elements are the variance of the noise in the measurement of the samples, and is a matrix of parameters of the model such that:(13) 
with a vector of parameters. The term acts as a regularization term that prevents overfitting and helps make the covariance matrix nonsingular.
The parameters that need to be learned are the elements of the vectors as well as the hyperparameters of the kernels .
In order to find the best hyperparameters of the model, the loglikelihood of the model on the data needs to be optimized. The solution of this nonconvex optimization problem is a subject of active research. Some software libraries rely on modelfree blackbox optimization techniques to solve it while others rely on gradientbased optimization techniques with multistart (to circumvent the nonconvexity), as is the case in the GPy [gpy] on which we rely in our work.
4.3. Optimization phase
Once the model of is either built (at the first iteration of ) or updated (at subsequent iterations), the optimization phase consists of finding the optima of an acquisition function over the model (such as EI), relatively to every task, through a modelfree blackbox optimization algorithm (such as Particle Swarm Optimization (PSO)). The resulting solution found for every task is then evaluated through an expensive run of the application to be tuned and then reinjected into the model.
The model can be queried for a prediction relative to a task and for an (unexplored) input , using not only the data of the task itself but also the data of the other tasks, through the formulation:
(14) 
A fundamental property of GPbased models (such as LCM) is the ability to estimate the confidence in the predictions alongside the predictions themselves. The confidence can be expressed as:
(15) 
There are several benefits to using a multioutput framework during the optimization phase.

The independence of the different optimizations (of the different tasks) allows for a high degree of parallelism, that can easily be exploited. This comes in addition to the parallelism available within the EGO algorithm itself.

The optimum of every task should be close to the optimum of related (close) tasks. Thus, the exploration around the optimum of one task will benefit the neighboring tasks as new data in the vicinity of their own optimum is gathered.

The effect of taking into account other tasks while predicting a given task increases the confidence in the prediction (decreases variance), which leads the to converge faster to the global optimum.
5. Transfer learning
In the field of machine learning, transfer learning consists of using the knowledge of one (or several) task(s) to improve the learning accuracy and / or speed of another task. Given the knowledge gathered and model built on previously autotuned representative tasks, we would like to be able to predict an optimal, or at least good enough, parameter configuration on a new unknown task. To this end, we propose two novel methods, the second one building on the first one. By combining in a new way some basic building blocks such as sampling techniques and GPs, they are able to target the needs of the multioutput autotuning framework that we develop in this paper.
The first method, that we denote as (short for Transfer Learning for Autotuning 1) and describe in Section 5.1, applies if an approximation to the optimum parameter configuration of a new task is considered enough, or, if a specific tuning for that task cannot be afforded.
The second method, that we denote as and describe in Section 5.2, applies if a high(er) quality of the optimal parameter configuration of a new task is desired and a quick additional tuning step can be afforded.
Let us note that both and benefit from each other. While relies on the approximation provided by to speed up its tuning, also benefits from the fact that updates the model described in Section 4.2 every time it is applied on a new task.
5.1. : Modeling the optima
The goal of the method that we propose is to create a model that can predict the optimal configuration corresponding to an unexplored task without having to tune it.
Given the tuning of the set of tasks , let us define the set of corresponding optima as
(16) 
An optimum parameter configuration is composed of as many parameters as the dimension of the input space , i.e. . Consequently, the solution that we propose is to create separate and independent Gaussian Processes to model every component of the optimal solutions separately. Such a Gaussian Process model is described in Section 4.2.1. Relating to , predicts and predicts , for example. The set represents the input data of every one of these Gaussian Processes. It is important to notice that, while the model in Section 4.2 has its space of inputs and space of outputs , every one of the Gaussian Processes created here has its space of inputs and its space of outputs in one of the dimensions of .
For any unexplored task , a prediction of its optimal parameter configuration is then given by
(17) 
Moreover, the confidence in the prediction of the models can serve as an indicator to the user of the autotuning on whether an additional tuning step (as described in the next section) is needed or not.
An alternative solution could have been to define a single multioutput Gaussian Process to model all the components simultaneously. However, no a priori hypothesis can be made in the general case on the correlations between the different components characterising the optima.
5.2. : Adding a new task
The goal of the method that we propose is to speed up the tuning of an unknown task by leveraging the knowledge of previously tuned tasks. It is composed of two steps.
5.2.1. Step 1: Centered sampling
The first step consists in predicting the optimum of the new task (through
) in order to focus the attention of the tuning to its surrounding area. Instead of sampling the whole search space of the new task with a general space filling sampling, a denser sampling is achieved in the area at the vicinity of the predicted optimum, while further away regions of the search space are less sampled. Specifically, given a budget of a certain number of initial samples, a normal distribution is used to map these samples on the search space. The center of the distribution is the predicted optimum. The standard deviation of the distribution is chosen as the diameter of the search space. Generated samples that are outside of the search space or violating constraints on the parameters are simply discarded and replaced by valid newly generated samples.
Let us note that, in the literature, the proposed sampling method is usually avoided. If used in conjunction with a regular Gaussian Process to model the objective function of the new task without any additional knowledge from other tasks, the EGO algorithm applied on this Gaussian Process will tend to favor the exploration of the undersampled areas, characterized by a high(er) prediction uncertainty. However, this would be unnecessary in our case as the model used in the second step helps circumvent this situation as the knowledge on the previously autotuned tasks informs the EGO algorithm of the low probability of finding an optimum in the low sampled regions of the space. Indeed, given the hypothesis of smoothness and continuity of the objective function between neighboring tasks, there should be no need to explore regions of the input space outside of the region neighboring the approximation to the optima.
A future improvement of the sampling is to chose a different variance of the distribution per dimension of the space. The lengthscale hyperparameter ( terms in Equation (5)) of the latent functions of the model of relative to every dimension of the search space can be used to figure out a relevant variance of the sampling distribution for that dimension. This improvement is simply inspired by Automatic Relevance Determination (ARD) methods [ard] where a relatively high value of the lengthscale hyperparameter relative to one dimension is an indicator of a small variance in that dimension, and viceversa.
5.2.2. Step 2: Extended model
The second step consists of extending the model of Section 4.2 to account for the new task to be tuned. Specifically, given the formulation of the model of a given task described in Equation (7), the matrices need to be extended by one more row and one more column. Given the formulation in Equation (5), this boils down to extending the vectors by a single additional scalar (hyperparameter), for all . This modified model has two computational advantages.
Firstly, instead of relearning all of the hyperparameters of the model, only the additional hyperparameters need to be learned, the other ones remaining unchanged. Indeed, the elements of the vectors describe the relation between the latent functions and the task models . Adding a new task should not affect this relation (too much) on the previously tuned tasks, and, even if it does, these tasks being already tuned, no learning is needed on them, and the transfer of the learning is meant to be only from them to the new task.
Secondly, the costliest part of the training and building of the model of both and is the inversion of the covariance matrix, as is the case with any Gaussian Processbased method. While a full factorization of the whole covariance matrix is needed in , efficient update / downdate techniques [Davis09] can be applied in . These techniques consist of reusing the factors of the previously factorized part of the covariance matrix and only inverting the newly added rows and columns of this matrix. In the case of , these new rows and columns correspond to the Kronecker product of the new rows and columns of the matrices with the covariance matrices generated by the kernels on the samples of the new task. The use of update / downdate techniques reduces the model building and training dramatically as they reduce the complexity of such operations from to .
5.3. with restarting
The guiding idea behind is that an expensive model is to be built once only, in on an initial set of (training) tasks, while every additional unencountered task can be tuned at a much reduced cost, by taking advantage of it and by only updating it rather than by recomputing it from scratch. However, this gain in model training cost comes at the expense of a potential loss in model accuracy.
While only adds and optimize new hyperparameters in the vectors, it does not modify the existing hyperparameters in these vectors, nor does it change the hyperparameters of the kernels related to the latent functions , nor does it increase the number of latent functions . In this regard, one can view the optimization strategy in as being analogous to stochastic gradient decent while that of would be analogous to gradient decent.
However, while it is a reasonable approximation to keep the old hyperparameters of the model fixed when adding few new tasks (as in ), it might become advantageous to update these hyperparameters if several new tasks are added, or even necessary to rebuild a new model from scratch, with a larger number of latent functions , if too many new tasks are added (as in ).
The choice of whether to continue relying on when a yet another new task is to be tuned rather than restarting , using instead, is highly dependant on: the application to be tuned; the tuning quality sought for; and the tuning cost a user is willing to pay. Specifically, two extreme regimes can be identified in autotuning. In the first extreme regime, such as tuning exascale application, every run on the objective function is very costly. Hence, the amount of data available to the tuner is small and every data is valuable. In such a context, the cost of optimizing the model is marginal. It is thus better for a user to restart every time a new task is to be tuned. In the second regime, every run of the application is cheap. Hence, the amount if data available is considerable. The cubic complexity of the inversion of the covariance matrix in the model would be a bottleneck to the tuning. It is then better for a user to rely on without any restart. For a given application, the need to restart is then depending on the cost of relative to the cost of a run of the application.
6. Experimental results
This section presents the experimental results of the methods described in this paper together with their comparison with stateoftheart autotuning methods. The autotuning of the QR factorization [demmel2012communication] routine in the ScaLAPACK library [slug] is explored as a guiding example.
We use the “Edison” machine at NERSC.^{1}^{1}1http://www.nersc.gov/users/computationalsystems/edison/ Each node of Edison contains dualsocket 12core Intel Ivy Bridge processors. The code was compiled with the Intel Fortran compiler version 18.0.1 and linked with craylibsci 18.03.1 for BLAS and ScaLAPACK operations.
The implementation of the algorithms proposed in this paper is done using the Python language. Also, the code relies on the GPy library [gpy] from University of Sheffield, which implements a general framework for Gaussian Processes.
Finally, to cope with runtime instabilities and noise in the measurements, all the runs of the application were performed times, and the minimal runtime was selected.
6.1. Analytical model
Although the methods in this paper do not rely on any analytical model but only rely on machine learning models, and, given the availability of such analytical model in the case of the QR factorization, we show here such models, although bringing insights about the behaviour of the application sought to be tuned, are not precise enough to be reliable models for autotuning, which highlights the importance of the use of blackbox machine learning models instead.
Detailed analytical models of QR factorizations runtime have been derived in the literature, some of which target distributedmemory algorithms. In this section, we show that such models are not accurate enough to predict performance.
In our study, we used three models. All of them estimate the runtime as a linear combination of factors and model parameters as follows:
(18) 
The factors represent the number of floating point operations , the number of messages and the volume of data transferred . The model parameters relative to each factor are the time per floating point operation , the time per message and the time per data item communicated
. They are hardware dependent and need to be estimated through a linear regression procedure.
The models differ in the varying levels of complexity and accuracy of the estimates of the factors they rely on. The first (basic) model is from [slug]. It relies only on the input matrix size and number of processes to compute the factors. The second (medium) model is formulated in [demmel2012communication]. It separates the effect of the number of divisions from the rest of the floating point operations. Its estimate of the factors is given by:
(19)  
where represents the block size (square blocks with ). The third (advanced) model given by [Laura] consists of using a MATLAB simulation that takes into account more detailed aspects of the factorization algorithm (particularly the [boundary conditions / edge cases]) in order to compute an accurate count of the factors in Equation (19).
In order to estimate the model parameters, we ran the QR factorization on a set of square matrices of size , and , and chose the values of the parameters following a semiexhaustive (grid) search. Given that the analytical models that we found in the literature consider only square blocks (), only experiments involving square blocks (between and instances) are used for the fitting, the remainder of the data (between and instances) is kept for model validation. We used the machine at TAMU for the experiments, which is composed of dual socket nodes, each socket containing 10 Intel Xeon E52670 v2 (Ivy BridgeEP) cores (2.5 GHz clock speed) and 64 Gb of memory. Updated version of Intel compiler and Intel MKL used for the experiments. ^{2}^{2}2https://hprc.tamu.edu/wiki/Ada
The average relative error is about for the basic model, and for the medium model when notseparating and when separating number of divisions from number of flops, respectively, and, for the advanced model. One reason for the large relative error is that these models make the assumption that most cores behave similarly and that the input matrices are large enough to run the level 3 BLAS routines at full speed. Although the predictions are not satisfying, our experiments show that there is a good correlation between experimental data and estimated data (of about ) for most experiments. Therefore the analytical models can still be used for autotuning purposes to show the trend of the objective function to optimize.
Moreover, we observed that a good set of parameters does not usually correspond to square blocks, although it is usually possible to find a set of parameters involving square blocks whose corresponding runtime is within of the runtime attained by the best parameters.
Furthermore, we observed that there are some patterns in the data and also that the accuracy of the analytical models can vary depending on the region of the search space. Further investigation is needed however to understand these effects.
6.2. Multitask learning autotuning
In this section, we compare three autotuning approaches: OpenTuner [opentuner], HpBandSter [bohb] and .
In order to be able to compare these methods, given that requires several tasks while the others only need one, we generated tasks following an LHS in the task space. The performance of a given tuning method on a given task is defined as the lowest computation time of the task yielded by the optimal (best) combination of parameters found by the tuning method. The stopping criteria of each method on every task is defined as a quota of evaluations of the objective function (runs of the application) that they are allowed to perform in order to find the optimum. In the case of , among these samples, the first ones are dedicated to the first phase (sampling phase) while the next ones are dedicated to the second and third phases of the method. This choice of the number of initial samples to be 3 times the number of parameters is a common practice in the Bayesian optimization literature. In general, it can affect the runtime and / or accuracy of the method. Also, OpenTuner is not able natively to take into account for constraints on the parameters. Similarly, HpBandSter can handle simple constraints but more advanced constraints cannot be handled, for example, the execution of a routine to check the validity of a combination of parameters. Thus, we reformulated the parameter space of the application in such a way as all the parameters configurations generated out of this modified space yield valid configurations in the initial parameter space.
Figure 1 presents comparative results between , OpenTuner and HpBandSter. Every point corresponds to a different task. Moreover, the sizes of the matrices considered is between and . Furthermore, compute nodes are used with cores per node. Any valid combination of MPI processes and threads is allowed. The Xaxis represents the ratio of best computation time of a task with OpenTuner (blue) and HpBandSter (orange) over that with . It compares the quality of the tuning of the different methods. A point at the right of the vertical line means that finds a better solution than the others, and viceversa. The Yaxis represents the ratio of the sum of the runtimes of all the samples generated by OpenTuner and HpBandSter over that with (excluding the time spent within the search algorithm itself). It compares the cost (in terms of application evaluations) of the different methods. A point above the horizontal line means that OpenTuner or HpBandSter are more expensive that , and viceversa.
Two overall conclusions can be drawn from this figure: () leads to better solutions than OpenTuner and HpBandSter as it leads to better application runtimes in () and () cases out of cases, respectively; () costs less in terms of total application runtime than OpenTuner in cases and HpBandSter in cases out of cases, respectively.
The average (up to ) improvement of the application runtime using compared to OpenTuner and HpBandster as observed in the Figure comes at the price of an increased tuning time. Indeed, although the cost in terms of total runtime of the application on the supercomputer is similar to that of the other tuners, the algorithmic complexity of is , where is the number of tasks and the number of samples per task. In contrast, the other tuners are executed times with a linear to quadratic algorithmic complexity in . Indeed, while OpenTuner and HpBandSter require on the order of to seconds to make a decision on the next parameter configuration to evaluate, it takes in the order of seconds for the model to be built (for the experiment of Figure 1) and an additional minute per task to find the best next parameter configuration. Since OpenTuner and HpBandSter tune every task separately, their total tuning cost is simply their above tuning cost per task times the number of tasks to be tuned. In contrast, the cost of building the model increasing cubically with the total number of runs of the application.
Fortunately, several parts of are amenable to parallelization. The costly inversion of the covariance matrix in the model can be achieved in parallel, either using the Cholesky factorization in a multithread LAPACK library, or using the distributedmemory Cholesky factorization in the ScaLAPACK library. Moreover, a second level of parallelism can be taken advantage of for the optimization of the hyperparameters of the model, either by using a parallel modelfree blackbox optimization technique, or by executing the different restarts of a gradientdecent based method in parallel. Furthermore, once the is built, the search for the best next parameter configuration can be done in parallel on the different tasks. Thus, when tuning exascale applications, the availability of large supercomputers is leveraged, hence alleviating the algorithmic complexity of .
Additionally, algorithmic improvements can be applied to reduce the cost of , potentially at a reduced model quality. Firstly, the number of latent functions is a parameter of the method that inversely impacts the quality and speed of learning of the model. Low values of lead to faster model learning but decrease the quality of the model. For instance, in our previous experiment, we have chosen the value of to be . Secondly, sparse approximations exist in the literature that approximate the covariance matrix (at the heart of the method) by a lowrank approximation which is much cheaper to deal with. Again, the lower the rank of the approximation, the worse the quality of the model. In the GPy package that we rely in the implementation of , the rank is influenced by the choice of a set of inducing points. In our previous experiment, we have chosen the number of inducing points as .
Hence, a balance can and must be found in the choice of and the rank / inducing points, which would depend on the tuning time that can be afforded together with the relative cost of a run of the application.
6.3. Transfer learning autotuning
Once the method has been applied on a set of tasks, leading to the construction of a performance model of the application, and are able to tune new tasks by transferring the knowledge acquired in to speed up and improve the accuracy of the tuning process.
This section compares the performance of OpenTuner, and . These three methods leading to better application runtimes than HpBandSter on all the cases tested in this experiment, its results are omitted here. new tasks are generated randomly in the task space. Similarly to the experiment in the previous section, OpenTuner and are given a quota of evaluations of the objective function, not needing any. Also, points are dedicated to the first phase of and others dedicated to the next phases. Moreover, the range of size of matrices and hardware configuration is the same as the ones used for in Section 6.2.
Figure 2 compares the different methods. The xaxis represents the different tasks. The yaxis represents the ratio of the best runtime for OpenTuner (blue bars) and (red bars) over that of . The horizontal line at represents the normalized runtime of . The tasks are ordered by increasing ratio of best runtime for OpenTuner over .
In this figure, the performance metric or success metric that we use to compare different tuning methods is the best runtime of the application found by a method on a given problem. From this figure, we can see that is competitive with OpenTuner as they outperform each other in of the cases. As does not perform any evaluation of the objective function, but guesses the optimal parameter configurations through a model prediction, this shows that the transfer of the learning from the model of is successful. Moreover, outperforms in cases, and performs similarly in the other cases. This result is not surprising as builds on top of , and incorporates the result of in its samplings. Finally, outperforms OpenTuner in cases out of and is outperformed in cases out of . is thus a viable competitor to stateoftheart autotuning methods.
6.4. SemiExhaustive search results
Figure 3 shows the performance of the routine of ScaLAPACK on the Edison computer, at NERSC, for two problem configurations, that is, a different , , . The computer being the same, . Figure 2(a) corresponds to and (sharedmemory setting), while Figure 2(b) corresponds to and (distributedmemory setting). On both its subfigures, the X and Y axes correspond to and , respectively. The Z axis represents the computation time. Given that the remaining parameters cannot be represented on a 3D graphic, every layer in the lasagnalike graphic correspond to a different combination of all of the other parameters (, , and ). Also, in every layer, the lighter the color, the higher the value of the objective function, while the darkest colors correspond to the regions of interest where the optima lies. The global optima correspond to the darkest region(s) at the lowest layer(s). The best parameters for each problem are described in Table 2:
Problem  Best parameters  Time (s)  

m  n  nodes  cores  mb  nb  nth  nproc  p  q  
2000  2000  1  24  4  8  1  24  2  12  0.20 
10000  10000  128  24  1  26  1  3072  16  192  1.61 
From these results, several conclusions can be drawn. Firstly, we can see that regions of interest exist, where the optima should be, and can be discovered by the optimization methods. Roughly speaking, they are located, in the first case, at small values of , while in the second case, at small values of both and . Secondly, the different lasagna layers are well separated in the first case while somehow interleaved in the second case. This last scenario makes it harder for an autotuner to figure out the optimal parameter configuration. Finally, in the second scenario, the regions where the optimum seems to be have a great amount of variability or noise in them, while noninteresting regions undergo much less variability. This can be explained by the fact that, small values of and lead to several small messages to be transferred over the network while other values of these parameters lead to fewer but larger messages to be transmitted. The first case is much more prone to be affected by network noise than the second case. Indeed, when large scale applications running on supercomputers, the applications share the computer network with each other. This phenomenon is a serious challenge for autotuning at the exascale.
If we consider that the best runtime found by means of grid search is close to optimal, it is possible to compare the other methods in terms of absolute performance (instead of relative). Table 3 offers such a comparison between grid search, OpenTuner, TLA1 and TLA2 on the tuning of a single task (matrix of size 500by500). Budget represents the number of runs of the application allowed (or necessary) for a specific method. mb, nb, nth, nproc, p and q represent the values of the parameters for the optimal configuration found by a given method, while the column Best Time shows the corresponding best runtime. OpenTuner finds the least good solution (compared to the other methods) when given a budget of and a good solution when given a budget of (impractical in an exascale setting). An approach such as that of OpenTuner is thus not adequate for tuning exascale applications, where the budget is limited. However, TLA1 (requiring 0 runs) and TLA2 (with a budget of ) find competitive solutions, even with such a low budget, making them more suitable candidates for tuning exascale applications.
Method  Budget  Time (s)  Best parameters  

mb  nb  nth  nproc  p  q  
OpenTuner  100  7.00e3  470  30  3  8  1  8 
OpenTuner  1000  5.30e3  88  12  1  24  1  24 
TLA1  0  5.73e3  244  11  1  24  1  24 
TLA2  100  5.67e3  296  11  1  24  1  24 
Grid search  8192  5.15e3  500  8  1  24  1  24 
7. Conclusion
This work presents the use of multitask and transfer learning for the purpose of autotuning. The core ideas are: ) multitask learning: instead of tuning an application on different input problems independently one from the other, the knowledge of all tasks can help tune every task faster, transfer learning: the knowledge accumulated from tuning previous tasks can be exploited to speed up the tuning of the new tasks.
The methods in this paper are intended to be used as follows. The user of an application would tune this application on an initial set of representative problems. This expensive tuning phase would be done once, through , in order to create a model of the objective function to optimize. Then, every time the application is to be tuned on a new task, either or would be used at a much reduced cost, the choice depending on the accuracy of the solution required or the amount of tuning time and computational resources available. If too many new tasks have been tuned through , the user might want to rebuild a new performance model using from scratch, as a large amount of data is available that was not during the initial build of the model.
The experimental results show that outperforms OpenTuner on of the cases they are compared on, with an application runtime improvement of up to . Moreover, the runfree and the lowcost methods find configurations whose runtimes outperform or are very competitive with those found by OpenTuner and HpBandSter. Furthermore, when the budget for autotuning (number of allowed runs of the application) is low, for instance, due to a high cost of the target application (e.g. exascale applications), our methods are more suitable than the other stateoftheart methods, which are not necessarily designed for this kind of use case.
This work presents a new view to autotuning that opens the doors to new challenges and research. Future directions include the use of: more complex models (e.g. deep Gaussian processes) to better take into account the nonstationarity, discontinuity and heavy noise which are characteristic behaviours of the objective functions measured on exascale applications (See Figure 2(b)), different models (e.g. Gaussian processes with random forests) to better handle the case of categorical parameters, and model approximations (e.g. sparse Gaussian processes) in order to scale the optimization algorithms themselves.