# Fast classical simulation of evidence for the utility of quantum computing
before fault tolerance

Tomislav Begušić  Garnet Kin-Lic Chan<EMAIL_ADDRESS>Division of Chemistry
and Chemical Engineering, California Institute of Technology, Pasadena,
California 91125, USA

###### Abstract

We show that a classical algorithm based on sparse Pauli dynamics can
efficiently simulate quantum circuits studied in a recent experiment on 127
qubits of IBM’s Eagle processor [Nature 618, 500 (2023)]. Our classical
simulations on a single core of a laptop are orders of magnitude faster than
the reported walltime of the quantum simulations, as well as faster than the
estimated quantum hardware runtime without classical processing, and are in
good agreement with the zero-noise extrapolated experimental results.

Clifford-based methods [1, 2, 3, 4, 5] are a promising approach for the
classical simulation of quantum circuits. The complexity of these approaches
is primarily governed by the number of non-Clifford gates because they use the
property that Clifford gates, which include all Pauli rotation gates with
angles $\theta=k\pi/2$ (for $k$ integer), can be simulated efficiently in
polynomial time. In contrast to other common techniques, such as tensor
network methods [6, 7, 8], Clifford-based approaches are not limited by the
amount of entanglement generated by the quantum circuit.

A recent experiment presented evidence for the utility of quantum computing
before fault tolerance by comparing results from the zero-noise extrapolated
[9] quantum simulation of the kicked Ising model on up to 127 qubits to those
from classical matrix product state (MPS) and isometric tensor network state
(isoTNS) simulations [10]. Ref. [10] showed that these tensor network
simulations failed to produce sufficiently accurate results for the set of
considered circuits, even in cases when all the rotation gates were Clifford
or near-Clifford. At the same time, the number of non-Clifford gates in the
experiment prohibit the use of exact exponentially-scaling Clifford approaches
[1, 2, 3, 4]. Here we apply a variant of a recently introduced Clifford
perturbation theory [5], which we term sparse Pauli dynamics (SPD), to the
simulation of circuits of Ref. [10]. The technique allows us to truncate the
exponential growing set of Pauli operators generated by the non-Clifford
gates, yielding an approximation algorithm that takes advantage of near-
Clifford structure.

Figure 1: Expectation values of different observables considered in Ref. [10].
Exact results [for panels (a)–(c)], quantum simulation experiment (“Quantum”),
and MPS results are reproduced from Figs. (3a)–(4b) of Ref. [10]. Following
that reference, the magnetization $M_{z}=\sum_{q}Z_{q}/n$ (a), weight-10 (b),
and weight-17 (c) observables are simulated for a quantum circuit consisting
of five Trotter steps (amounting to five layers of non-Clifford gates),
another weight-17 observable (d) is computed for a quantum circuit consisting
of five Trotter steps and an additional $R_{X}$ layer (six non-Clifford layers
in total), and $\langle Z_{62}\rangle$ in (e) is simulated for a circuit
composed of twenty Trotter steps (twenty non-Clifford layers). For visual
simplicity, values in (c) and (d) were negated. The SPD results were computed
at $\theta_{h}=k\pi/32$ for integer $k$ between 0 and 16 (“+” markers).

Results. Following Ref. [10], we consider a quantum circuit composed of
alternating layers $R_{ZZ}=\prod_{\langle
i,j\rangle}\exp(i\frac{\pi}{2}Z_{i}Z_{j})$ and
$R_{X}(\theta_{h})=\prod_{i}\exp(-i\frac{\theta_{h}}{2}X_{i})$ of Pauli
rotation gates, i.e.,

$U=[R_{ZZ}R_{X}(\theta_{h})]^{t},$ (1)

where $\langle i,j\rangle$ denotes that $i$ and $j$ qubits are neighbors in
the heavy-hexagon lattice of IBM’s Eagle processor [10] and $t$ is the number
of Trotter steps. Since the entangling gates in $R_{ZZ}$ are Clifford, we
first transform the circuit to consist only of non-Clifford Pauli rotation
gates $R_{\tilde{X}^{(m)}}$, where
$\tilde{X}_{i}^{(m)}=(R_{ZZ}^{\dagger})^{m}X_{i}R_{ZZ}^{m}$. We can further
perform an angle transformation of non-Clifford rotation gates
$R_{X}(\theta_{h})=R_{X}(\tilde{\theta}_{h})R_{X}(\pi/2)$. We apply this
transformation in the simulation of weight-10 and weight-17 observables (see
below), which under this transformation convert into single-qubit Pauli
operators. To handle the remaining non-Clifford gates, we then perform the
sparse Pauli dynamics, which works in the Heisenberg picture of the observable
in the Pauli representation (see Method). Its accuracy is controlled by a
truncation order $K$, which limits the size of the Pauli basis.

Figure 1 compares our classical simulation results with those obtained in Ref.
[10] by a numerically exact method [for panels (a)–(c)], the quantum
simulation experiments, and the MPS approach. The maximum truncation orders
were $K=10$ (a, b, e) and $K=6$ (c, d). With these parameters, each point in
Fig. 1 takes between one and two minutes on a single processor of a laptop
computer, which is two orders of magnitude faster than the reported quantum
wall-clock run time for Figs. 1d, e (4 h and 9.5 h, respectively) and faster
than the hypothesized run time of IBM’s Eagle processor without any classical
processing steps (estimated at about 5 min for Fig. 1e). Although the chosen
truncation orders do not yield numerically exact results (which can be easily
verified by the small deviations in panels (a)–(c)) and more detailed studies
of convergence should be carried out, the results are comparable to (and
generally well within) the estimated errors of the quantum simulation
experiment and much more accurate than the MPS results of Ref. [10].

Overall, our results demonstrate that our simple sparse Pauli dynamics
approach can faithfully simulate the quantum experiments recently performed on
the new Eagle processor of IBM. Two days before this preprint was posted, a
preprint was posted by Tindall et al. [11] that showed that an alternative
tensor network method that takes advantage of the qubit connectivity of the
device also outperforms the MPS or isoTNS [12] methods employed in Ref. [10]
and can also simulate the quantum experiments. In contrast to that work, our
method is not specially constructed for a given qubit connectivity map, nor
does it rely on low entanglement. That both classical approaches succeed
illustrates the rich landscape of approximate classical algorithms that remain
to be explored. Finally, we believe that the method we describe here holds
promise not only for quantum circuit simulations, but also for more general
simulation problems in quantum dynamics.

Method. Consider the unitary composed of Pauli rotation gates after the
Clifford transformation described above:
$U=\prod_{i}U_{j}(\theta_{j})=\prod_{j}\exp(-i\theta_{j}P_{j}/2)$, defined by
the rotation angles $\theta_{j}$ and Hermitian Pauli operators $P_{j}$. (Note
that because the Clifford gates have been applied, the Paulis $P_{j}$ are not
those of the original circuit, and all the rotations are non-Clifford).

We compute the expectation value of a Pauli observable

$\langle O\rangle=\langle 0^{\otimes n}|U^{\dagger}OU|0^{\otimes n}\rangle$
(2)

by performing Heisenberg evolution of the observable under Pauli rotation
gates, using

$U_{j}(\theta_{j})^{\dagger}\sigma
U_{j}(\theta_{j})=\cos(\theta_{j})\sigma+i\sin(\theta_{j})P_{j}\sigma$ (3)

whenever $\\{\sigma,P_{j}\\}=0$ and $U_{j}(\theta_{j})^{\dagger}\sigma
U_{j}(\theta_{j})=\sigma$ otherwise. For Clifford gates, Eq. (3) simply
transforms a Pauli operator $\sigma$ into another Pauli operator, whereas for
non-Clifford gates it results in a sum of two Pauli operators, leading to an
exponential computational cost. This sum can be written as

$\displaystyle\langle O\rangle$ $\displaystyle=\sum_{k=0}^{N}E^{(k)},$ (4)
$\displaystyle E^{(k)}$ $\displaystyle=i^{k}\sum_{1\leq
j_{1}<j_{2}<\cdots<j_{k}\leq
N}c_{j}\sin(\theta_{j_{1}})\cdots\sin(\theta_{j_{k}})\langle 0^{\otimes
n}|P_{j_{k}}\cdots P_{j_{1}}O|0^{\otimes n}\rangle,$ (5)

where $N$ is the number of gates, $j_{1},\dots,j_{k}$ are selected so that
$\\{P_{j_{m}},P_{j_{m-1}}\cdots P_{j_{1}}O\\}=0$, and $c_{j}$ is a product of
cosines of $\theta_{l\notin\\{j_{1},\dots,j_{k}\\}}$ [5]. To truncate the
exponentially growing sum of Pauli operators, we recently introduced a
perturbation-inspired expansion [5], which assumes sufficiently small angles
$\theta_{j}$ [i.e., small $\sin(\theta_{j_{m}})$] that suppress high-order
terms $E^{(k)}$. This assumption is supported by the fact that the rotation
gates can be rewritten as
$U_{j}(\theta_{j})=U_{j}(\tilde{\theta}_{j})U_{j}(k\pi/2)$ to ensure that
$\tilde{\theta}_{j}\in[-\pi/4,\pi/4]$, while $U_{j}(k\pi/2)$ are Clifford and
do not directly contribute to the computational cost.

In contrast to our previous work, here we do not directly truncate the series
(4) up to a specified order $K$ but use the same order concept to adaptively
increase the Hilbert (operator) subspace on which the approximate Heisenberg-
evolved operator is supported. Because the Heisenberg observable is supported
in a sparse subspace of the Pauli basis, we refer to this scheme as sparse
Pauli dynamics. For related concepts, see also [13, 14]. Specifically, let

$O_{j}=\sum_{\sigma\in\mathcal{P}_{j}}a_{\sigma,j}\sigma$ (6)

be the observable evolved with $j$ gates, $\mathcal{P}_{j}$ a set of Pauli
operators at step $j$, and $a_{\sigma,j}$ complex coefficients. In addition,
we assign an order $k_{\sigma}$ to each Pauli operator $\sigma$, defined as
the smallest number $k$ for which $\sigma=P_{j_{k}}\cdots P_{j_{1}}O$ (note
that $P_{j_{i}}$ is not in general a single qubit Pauli, so the weight of the
operator string is in general $>k+1$). The rotation gate
$U_{j+1}(\theta_{j+1})$ will commute with a subset of Pauli operators in
$\mathcal{P}_{j}$ and will transform the rest according to Eq. (3). Then, the
coefficients are updated as

$a_{\sigma,j+1}=\cos^{\xi_{\sigma,P_{j+1}}}(\theta_{j+1})a_{\sigma,j}+i\xi_{\sigma,P_{j+1}}\sin(\theta_{j+1})a_{P_{j+1}\sigma,j},$
(7)

where $\xi_{\sigma,\sigma^{\prime}}$ is 1 if the operators anticommute and 0
otherwise. For compactness, in Eq. (7) it is understood that $a_{\sigma,j}=0$
if $\sigma\notin\mathcal{P}_{j}$, even though in practice this zero
coefficient is not stored. The number of Pauli operators increases by one
whenever $\xi_{\sigma,P_{j+1}}=1$, $\sigma\notin\mathcal{P}_{j}$, and
$\sigma^{\prime}=P_{j+1}\sigma\in\mathcal{P}_{j}$, because then we must add
$\sigma$ in the set $\mathcal{P}_{j+1}$. To control the number of Pauli
operators from growing exponentially, we discard $\sigma$ with $k_{\sigma}>K$,
where $K$ is the maximum truncation order and a parameter controlling the
accuracy of the simulation. Finally, the expectation value is obtained as

$\langle 0^{\otimes n}|O_{N}|0^{\otimes
n}\rangle=\sum_{\sigma\in\mathcal{P}_{N}}a_{\sigma,N}\langle 0^{\otimes
n}|\sigma|0^{\otimes n}\rangle.$ (8)

A Python implementation using IBM’s Qiskit [15] was employed in producing the
results for Fig. 1.

Acknowledgments. TB and GKC were supported by the US Department of Energy,
Office of Science, Office of Advanced Scientific Computing Research and Office
of Basic Energy Sciences, Scientific Discovery through Advanced Computing
(SciDAC) program under Award Number DE-SC0022088. TB acknowledges financial
support from the Swiss National Science Foundation through the Postdoc
Mobility Fellowship (grant number P500PN-214214). GKC is a Simons Investigator
in Physics.

## References

  * Bravyi and Gosset [2016] S. Bravyi and D. Gosset, Phys. Rev. Lett. 116, 250501 (2016).
  * Bravyi  _et al._ [2019] S. Bravyi, D. Browne, P. Calpin, E. Campbell, D. Gosset, and M. Howard, Quantum 3, 181 (2019).
  * Bennink  _et al._ [2017] R. S. Bennink, E. M. Ferragut, T. S. Humble, J. A. Laska, J. J. Nutaro, M. G. Pleszkoch, and R. C. Pooser, Phys. Rev. A 95, 062337 (2017).
  * Qassim, Wallman, and Emerson [2019] H. Qassim, J. J. Wallman, and J. Emerson, Quantum 3, 170 (2019).
  * Begušić, Hejazi, and Chan [2023] T. Begušić, K. Hejazi, and G. K.-L. Chan, “Simulating quantum circuit expectation values by Clifford perturbation theory,” (2023), arXiv:2306.04797 [quant-ph] .
  * Gray and Kourtis [2021] J. Gray and S. Kourtis, Quantum 5, 410 (2021).
  * Huang  _et al._ [2021] C. Huang, F. Zhang, M. Newman, X. Ni, D. Ding, J. Cai, X. Gao, T. Wang, F. Wu, G. Zhang, _et al._ , Nature Computational Science 1, 578 (2021).
  * Gray and Chan [2022] J. Gray and G. K.-L. Chan, “Hyper-optimized compressed contraction of tensor networks with arbitrary geometry,” (2022), arXiv:2206.07044 [quant-ph] .
  * Temme, Bravyi, and Gambetta [2017] K. Temme, S. Bravyi, and J. M. Gambetta, Phys. Rev. Lett. 119, 180509 (2017).
  * Kim  _et al._ [2023] Y. Kim, A. Eddins, S. Anand, K. X. Wei, E. Van Den Berg, S. Rosenblatt, H. Nayfeh, Y. Wu, M. Zaletel, K. Temme, _et al._ , Nature 618, 500 (2023).
  * Tindall  _et al._ [2023] J. Tindall, M. Fishman, M. Stoudenmire, and D. Sels, “Efficient tensor network simulation of IBM’s kicked Ising experiment,” (2023), arXiv:2306.14887 [quant-ph] .
  * Zaletel and Pollmann [2020] M. P. Zaletel and F. Pollmann, Phys. Rev. Lett. 124, 037201 (2020).
  * Kuprov, Wagner-Rundell, and Hore [2007] I. Kuprov, N. Wagner-Rundell, and P. Hore, Journal of Magnetic Resonance 189, 241 (2007).
  * Rakovszky, von Keyserlingk, and Pollmann [2022] T. Rakovszky, C. von Keyserlingk, and F. Pollmann, Physical Review B 105, 075131 (2022).
  * Qiskit contributors [2023] Qiskit contributors, “Qiskit: An open-source framework for quantum computing,”  (2023).


# New approach to MPI program execution time prediction
††thanks: This work is supported by Russian Ministry of Science and Higher
Education, grant N 05.613.21.0088, unique ID RFMEFI61318X0088 and the National
Key R&D Program of China (2017YFE0123600)

Andrey Chupakhin  Lomonosov Moscow State University
Moscow, Russian Federation
<EMAIL_ADDRESS>Alexey Kolosov  Lomonosov Moscow State University
Moscow, Russian Federation
<EMAIL_ADDRESS>Ruslan Smeliansky  Applied Research Center for
Computer Networks (ARCCN)
Lomonosov Moscow State University
Moscow, Russian Federation
<EMAIL_ADDRESS>Vitaly Antonenko  Applied Research Center for
Computer Networks (ARCCN)
Lomonosov Moscow State University
Moscow, Russian Federation
<EMAIL_ADDRESS>Gleb Ishelev  Lomonosov Moscow State University
Moscow, Russian Federation
<EMAIL_ADDRESS>
###### Abstract

The problem of MPI programs execution time prediction on a certain set of
computer installations is considered. This problem emerge with orchestration
and provisioning a virtual infrastructure in a cloud computing environment
over a heterogeneous network of computer installations: supercomputers or
clusters of servers (e.g. mini data centers). One of the key criteria for the
effectiveness of the cloud computing environment is the time staying by the
program inside the environment. This time consists of the waiting time in the
queue and the execution time on the selected physical computer installation,
to which the computational resource of the virtual infrastructure is
dynamically mapped. One of the components of this problem is the estimation of
the MPI programs execution time on a certain set of computer installations.
This is necessary to determine a proper choice of order and place for program
execution. The article proposes two new approaches to the program execution
time prediction problem. The first one is based on computer installations
grouping based on the Pearson correlation coefficient. The second one is based
on vector representations of computer installations and MPI programs, so-
called embeddings. The embedding technique is actively used in recommendation
systems, such as for goods (Amazon), for articles (Arxiv.org), for videos
(YouTube, Netflix). The article shows how the embeddings technique helps to
predict the execution time of a MPI program on a certain set of computer
installations.

###### Index Terms:

Pearson correlation coefficient, matrix decomposition, embeddings, MPI,
execution time prediction, ensemble

##  I Introduction

The idea of building a virtual infrastructure in a cloud computing environment
over a heterogeneous network of computer installations is not new [1], [3],
[7], [8]. To do this several administrative authorities called federates,
allocate some their computational resources for the cloud computing
environment operation. Such unions, called federations, have already been
established [1], [4], [5], [6], some are only planning to be created [2], [9].
Each federate can consists of several computer installations with different
computer powers, such as high-performance computing cluster, data center or
supercomputer. Federate provides resources of its computing environment
111Further the computing environment will be understood as a heterogeneous
network of computer installations for creating a virtual infrastructure in
accordance with the access policy set by the federate administration. Further
in this article, the term program will be understood as MPI Program.

One of the main criteria for the effectiveness of such cloud computing
environment is the time spending by a program in this environment. This time
consists of the time spent by program in the queue (waiting time) and the
program execution time (execution time). This value depends on the resources
allocation algorithms of the cloud environment (CE) (mapping virtual computer
installations to physical ones) and the queue discipline, taking into account
the heterogeneity of physical computer installations.

Consider the “path” of a program from the entering the CE up to getting the
result of its execution:

  1. 1.

The user forms a set of input data: the text of the program, the program
arguments, program launch arguments (e.g. the number of MPI processes or the
requested resources (cpu, ram)), the input data for program, a special script
to compile program;

  2. 2.

The user send a set of input data to the CE through a single portal of the
federation. The input data from the portal come to the orchestrator of the CE,
which is responsible for implementing the virtual infrastructure and executing
the programs in it;

  3. 3.

The orchestrator responsible for federation resources allocation and has a
unique algorithm to do this. This algorithm determines the best order and the
best federate and physical computer installation inside federate for program
execution, in sense of the above mentioned effectiveness criterion. When
computer installation is selected this program is sent to there;

  4. 4.

Each federate has its own queue of programs. The arrived program is processed
by the local resources scheduler and will start running at a certain time
defined by the local scheduler;

  5. 5.

When the program execution is completed the results is returned to the
orchestrator which in its term sends it in the proper form to the users.

From the path description above, it can be seen that step 3 (meta scheduling)
and step 4 (local scheduling) can be optimized in the context of this article.
The focus of this article is step 3, namely: how to choose the most effective
computer installation that is a part of the federation computing environment
for the received program based on the minimum execution time criterion. The
article offers algorithms for program execution time prediction on a certain
set of computer installations based on the execution history of the program on
different computer installations.

The program execution time on the specific computer installation, as well as
the waiting time in the queue, can be predicted based on the history of its
launches on this computer installation [22], [23], [24], [25]. To do this,
various extrapolation algorithms can be used [26], [27], regression [10], or
more complex algorithms are used, such as the ensemble of decision trees
(Random Forest) [11]. The main disadvantage of these algorithms is that they
applicable only to a single certain computer installation. But the point of
the question of this article is a program execution time prediction on a
certain set of computer installations. These are the computer installations
whose characteristics meet the requirements of the virtual infrastructure. Of
course, the algorithms that were mentioned above can be used to get estimates
for programs execution time on the several computer installations. However,
this requires a history of running of this program on each computer
installation from the set. As a rule, this information isn’t available.

The logic for choice of computer installation based on the histories of
program execution can be described as follows. One of the well-known
algorithms [22], [23], [24], [25] is applied to program execution histories to
estimate program execution time on each of the computer installation from a
certain set of computer installations. With obtained estimates, one can either
developing scheduling for a group of programs, or follow a greedy strategy and
send each program to the computer installation that has the minimum execution
time.

As noted above, in order to allocate programs between several computer
installations, a history of program executions on each of the computer
installations is required.

As demonstrated in this paper this requirement can be relaxed: to predict the
program execution time on several computer installations, some running
histories of this program on them are sufficient. Moreover, it is not
necessary that each program launches on each computer installation. The
accuracy of the prediction depends on the number of program running histories
on computer installations from a certain set. To do this, the article proposes
a new approaches to estimate program execution time, which allows predicting
the program execution time on a certain set of computer installations. The
algorithms presented in this article for implementing these approaches allow
one to predict the execution time of a program on a certain computer
installation, even there is no execution history the program on it.

In the proposed approaches, the problem of the program execution time
prediction is reduced to filling empty entries in the matrix “Programs-
Computers” built for a given set of programs and a given set of computer
installations, in the entries which is an execution time a certain program
with certain sets of arguments on a specific computer installation. For more
information about the mathematical statement of this problem see Section 2.

In this article proposed two approaches to the program execution time
prediction problem on a certain set of computer installations. The first one
is based on computer installations grouping based on the Pearson correlation
coefficient [31]. It is shown that this method should be used in the case of a
”densely” (dense matrix) filled matrix ”Programs-Computers” (at least $95\%$
of entries are filled). The second one is based on decomposition of “Programs-
Computers” matrix into vector representations of computer installations and
programs, so-called embeddings. Here embedding is treated as a structured type
of unstructured data. For example, the form of a vector implies the structure
that is important here in contrast to hashing, which simply converts a
potentially infinite object into a finite one. In the paper it is shown how to
use embedding of program and embedding of computer installation to predict the
program execution time on specific computer installation. To calculate
embeddings, it’s proposed to use the technique [Matrix-Decomposition]
decomposition of the matrix ”Programs-Computers”. In Section 6 it is shown by
experimental studies that the second method works well with sparse matrices.
Histories of runnings of MPI and OpenMP benchmarks [33], [34] on a dozens of
different computer installations were used for experiments. This data is
described in Section 4.

It’s important to emphases that, as it will be shown below, the proposed
decomposition technique results embeddings of dimension 1, that allows to
place a total order relation on a set of computer installations and programs.

The rest part of the article is structured as follows. The 5th Section
contains a description of the developed algorithms for programs execution time
prediction: an algorithm based on computer installations grouping that works
well for dense matrices and an algorithm for obtaining vector representations
(embeddings) of computer installations and programs that works well for sparse
matrices. Section 7 describes the application of the developed algorithms to
solve the problem of ”minimization of the time spent by the program inside the
cloud environment”. In the last 8-th section, conclusions are presented and
goals for further research are formulated.

##  II Problem Description

We consider the case when almost nothing is known about the program and the
computer installation. The program source code and binary code are not
available, the architecture of the computer installation is unknown. Only the
program execution time, its arguments (see point 1 in Section 1) and the
numbers of used resources of computer installation are known.

Further, the program will be considered in the form of job which has the
following characteristics: the time when job is submitted to the queue, the
start time and end time, resource utilization (cpu, ram, disk), input data,
arguments of the program and arguments of the computing environment in which
this job is started (e.g., the number of MPI processes, the number of threads
for OpenMP, and the amount of resources needed to complete the job).

It is assumed that only the total amount of resources on the computer
installation is known: the number of nodes/CPUs/cores, the amount of RAM and
the storage space.

We will use the following notations to formulate the problem:

  1. 1.

$P_{i}$ – unique identifier of program;

  2. 2.

$\\{P_{1},P_{2},\ldots,P_{N}\\}=\\{P_{i}\\}$ – a set of $N$ unique identifiers
of programs; $Arg$ – arguments of program: program’s arguments (for the
executable file) and arguments of the computing environment where the program
is running (the number of MPI processes, the amount of requested computing
nodes /CPUs /cores);

  3. 3.

$\\{Arg_{1},Arg_{2},\ldots,Arg_{A_{i}}\\}=\\{Arg_{j}\\}$ – here each $Arg_{i}$
is the set of arguments of $i$-th program. A total amount of different sets of
arguments – $A_{i}$;

  4. 4.

$[(P_{i},Arg_{1}),(P_{i},Arg_{2}),\ldots,(P_{i},Arg_{A_{i}}),(P_{i},Arg_{1})]$
– a history of $P_{i}$ program running, arguments can be the same in some
pairs;

  5. 5.

$C_{i}$ – unique identifier of computer installation;

  6. 6.

$\\{C_{1},C_{2},\ldots,C_{M}\\}=\\{C_{i}\\}$ – a set of $M$ unique identifiers
of computer installations;

  7. 7.

Matrix $PC$ ($P$ – “program”, $C$ – “computer”) – the matrix where rows
correspond to history of running of programs, columns correspond to the
computer installations. This matrix is divided into groups of rows. There are
$N$ groups in total (by the number of the programs). The $i$-th group
corresponds to the history of running of the $i$-th program. The group
consists of $A_{i}$ rows, where each matrix entry display the execution time
of the program $P_{i}$ with the corresponding arguments $Arg_{i}$ on the
corresponding computer installation. The matrix entry is empty if the program
was never run on this computer installation.

TABLE I: Example of $PC$ matrix $P/C$ | $C_{1}$ | $C_{2}$ | $C_{3}$ | $C_{4}$
---|---|---|---|---
$P_{1},Arg_{1}$ | ? | ? | 90 | ?
… | … | … | … | …
$P_{1},Arg_{A_{1}}$ | 5 | 10 | 15 | 25
… | … | … | … | …
$P_{i},Arg_{1}$ | 25 | 35 | ? | 56
… | … | … | … | …
$P_{i},Arg_{A_{i}}$ | 60 | 75 | 96 | ?
… | … | … | … | …
$P_{N},Arg_{1}$ | ? | 34 | 67 | ?
… | … | … | … | …
$P_{N},Arg_{A_{1}}$ | 100 | 146 | 245 | 300

In the Table I an example of a $PC$ matrix is presented. Empty entries are
marked with a question mark. The $PC$ matrix shows the history of finished
program executions. It can be built from log files of the computer
installations – $C_{1},...,C_{M}$.

The terms of the notations introduced above the problem of program execution
time prediction on a set of computer installations can be formulated as
following:

  * •

Given

    * –

A set of $N$ unique identifiers of programs – $\\{P_{i}\\}$;

    * –

For each program $P_{i}$ a set with $A_{i}$ sets of arguments –
$\\{Arg_{j}\\}$, $|\\{Arg_{j}\\}|=A_{i}$;

    * –

A set of $M$ unique identifiers of computer installations – $\\{C_{i}\\}$;

    * –

A $PC$ matrix where entries contain the execution times of programs $P$ (with
the corresponding arguments – $Arg$) on computer installation $C$;

  * •

To find

    * –

Fill in empty entries of the $PC$ matrix, or in other words predict the
program execution times on the computer installations.

To evaluate the program execution time prediction, the prediction error is
calculated as follows:

$PredictionError(P_{i},Arg_{j})=\frac{|predict-target|}{target}$ (1)

where predict is predicted time of program $P_{i}$ with arguments $Arg_{j}$,
target is a true execution time of program $P_{i}$ with arguments $Arg_{j}$.

The total prediction error is calculated as an average of the errors
calculated using the equation (1) for all programs. This allows you to compare
algorithms among themselves.

##  III Related Works

The problem of the execution time prediction has been studied for a very long
time. For example, for real-time systems, the problem of estimation the worst-
case execution time (WCET) is still one of the main problems, since WCET
programs are used to build a schedule of tasks in real-time systems. In [13],
the execution time of specific applications on GPUs is predicted. In [14] a
model of the program and the computer are built in order to estimate the
program execution time on the computer. In [15] the execution time of programs
on Grid systems are predicted. For time prediction purposes are used
Analytical approaches [16], [17], [18], statistical approaches [19], [20],
[21] approaches based on historical data [22], [23], [24], [25], time series
prediction [26], [27], neural networks are also used [28], [29]. Execution
time can be predicted from test runs [30].

All algorithms for predicting program execution time mentioned above use the
history of program executions to make a prediction. Their main drawback is
that these algorithms only work when the history of program executions is
known on the computer installation, i.e. to predict the execution time on a
certain computer installation, you need more than one run of the program on
this computer installation because the history of executions have to consists
of at least two executions.

Besides above, in order to build a program model or make test runs, you need
the source code and executable files, and this information is often not
available. In addition, almost all articles that consider the prediction of
the program execution time often considered only for one computer
installation. However, in the point of question of this paper it has to be
able to estimate the program execution time on a set of computer installations
in that case when there is a minimum amount of information about computer
installations and programs, actually available.

##  IV Data for Prediction Model

Information about programs (benchmarks) executions on many computer
installations from the website spec.org was used for the experiments presented
in this article. This information was used to form the matrix “Programs-
Computers” $PC$, described in Section 2 above. Two data sets with the
execution results of programs on different computer installations are
presented in [33]:

  1. 1.

MPIL2007 – $12$ programs, $163$ computer installations;

  2. 2.

MPIM2007 – $13$ programs, $396$ computer installations.

It should be taken into account that the data on the website spec.org is
constantly updated. The data are used in the presented research dated April 8,
2020 12:13. The history of changes can be viewed in [35]. The execution
results of OpenMP benchmarks running [34] were also used as a data to form the
$PC$ matrix.

The problem formulated in Section 2 implies that the program can be run with
different arguments. However for the simplicity reason we will assume that
each program was run only with one set of arguments, i.e. one program
corresponds to one row in the $PC$ matrix. Note that this assumption does not
limit the generality of the proposed approaches.

A brief description of the data sets is provided in the Table II. Each entry
in [33] was considered as a separate computer installation.

TABLE II: Results of MPI and OpenMP programs launches on various computer installations | Number of | Number of | Benchmark |
---|---|---|---|---
Name | computer | programs | type | Link to data
| installations |  |  |
MPIL2007 | 163 | 12 | MPI | [33]
MPIM2007 | 396 | 13 | MPI | [33]
ACCEL_OMP | 25 | 15 | OpenMP | [34]

The used data was uploaded on
Github222https://github.com/andxeg/CloudCom_2020_Embeddings [36].

##  V Proposed Solutions

To solve the problem of the programs execution time prediction (see Section 2)
the following algorithm was developed (for brevity, only the algorithm scheme
is presented):

  1. 1.

Data preparation. Form a $PC$ matrix of dimension $N$ by $M$;

  2. 2.

Execution time prediction. Applying one of the algorithms described below:

    1. (a)

Ridge. Using ridge regression [10] for each row of the $PC$ matrix. In other
words the unknown execution time on some computer installation was predicted
from the known execution times of the particular program on all other computer
installations;

    2. (b)

Cliques. Grouping computer installations using the Pearson correlation
coefficient;

    3. (c)

Decomposition. Apply matrix decomposition [32] or tensor decomposition [46] of
the $PC$ matrix in order to fill in empty entries of the $PC$ matrix;

    4. (d)

Apply an ensemble of algorithms: Ridge \+ Cliques \+ Decomposition;

  3. 3.

Evaluation of prediction quality. The prediction of program execution times is
evaluated as follows:

    1. (a)

For each program a prediction error is calculated by the equation (1)
specified in Section 2;

    2. (b)

The error average is a general prediction error.

As it was mentioned above each program was run only with one set of arguments,
i.e. one program corresponds to one row in the $PC$ matrix. If the program was
executed on multiple sets of arguments, it’s mapped to multiple rows in the
$PC$ matrix. In this case, one can first use the history-based execution time
prediction algorithm [22], [23], [24], [25] to fill in the empty entries, and
then use the algorithm described above in 2a-2d. For more information about
the case when a single program corresponds to multiple rows in the $PC$
matrix, see Section 7.

###  V-A Computer installations grouping based on the Pearson correlation

The Pearson correlation coefficient [31] characterizes the presence of a
linear relationship between two sets of numbers. The columns of the $PC$
matrix are considered as such sets. By calculating the correlation between the
columns $C_{i}$ and $C_{j}$, one can detect a linear relationship between the
programs execution times on the corresponding computer installations. If the
Pearson correlation between columns $C_{i}$ and $C_{j}$ is close to $1$ in
modulus, then there is a linear relationship between the execution times of
programs on $C_{i}$ and $C_{j}$ computer installations. Therefore, a column
with program execution times on the $C_{i}$ computer installation can be
obtained by multiplying by a certain constant the column with program
execution times on the $C_{j}$ computer installation.

The grouping procedure for computer installations consists of the following
steps:

  1. 1.

Calculate Pearson correlation for each pair of columns of the matrix $PC$ –
total $\frac{M(M-1)}{2}$ pairs, where $M$ is the amount of computer
installations;

  2. 2.

Build the graph of computer installations:

    1. (a)

Each vertex represents one computer installation;

    2. (b)

An edge between the $C_{i}$ and $C_{j}$ computer installations exists if the
Pearson correlation between the corresponding columns in the $PC$ matrix is
modulo greater than some threshold. The value of threshold is an algorithm
parameter;

  3. 3.

Find all cliques [37] in the graph of computer installations. Cliques
searching is NP-complete problem. To solve this problem you can use the
algorithm from [38]. To speed up the experiments, a simplified cliques search
algorithm was implemented [36];

  4. 4.

Each clique is a group of computer installations with execution times of each
program are linear related, i.e. each pair of computer installations in such a
group has execution times related by a specific scalar coefficient. Pay
attention that the same computer installation can fall into several groups;

  5. 5.

The resulting cliques are groups of computer installations.

To search for cliques, the simplified algorithm was developed with complexity
no bigger than $N^{3}$ [36], which is significantly less than the complexity
of the algorithm for searching for all maximum size cliques [38] –
$3^{\frac{N}{3}}$. It finds at least one clique for each vertex in the graph.
The disadvantage of the proposed algorithm is that it doesn’t find all the
maximum cliques. However, if threshold for the value of Pearson correlation is
close to $1$, then further prediction algorithm described below is pretty good
even some vertexes in the cliques would be missed.

Based on computer installations grouping, the procedure for predicting program
execution time on the computer installation can be described as follows:

  1. 1.

Select $P$ – the program for execution time prediction;

  2. 2.

Select $C$ – the computer installation where $P$ program will be run;

  3. 3.

Determine the group that $C$ belongs to;

  4. 4.

If the computer installation doesn’t belong to any group of computer
installations, then make a prediction using the Ridge Regression (See Section
5, point 2.a). Regression is applied for the entire row in $PC$ matrix
correspond to the $P$ program;

  5. 5.

If the computer installation $C$ belongs to a certain group, then there is a
clique in the graph of computer installations, in which all the computer
installations, including $C$, are connected to each other. This means that
Pearson correlation is greater than some given threshold for each pair
(computer installation $C$; some other computer installation from the clique –
$C^{\prime}$). Therefore, it’s possible to calculate the coupling coefficient
between the programs execution times for two computer installations from a
pair;

  6. 6.

Then multiply the $P$ program execution time on the computer installation
$C^{\prime}$ on this coefficient. As result we get the estimate of $P$ program
execution time on the considered computer installation $C$;

  7. 7.

This procedure should be applied to each computer installation from the
clique, and then calculate the average of the execution times. The average
time is the estimation of the $P$ program execution time on the $C$;

  8. 8.

If it isn’t possible to calculate Pearson correlation (e.g, the considered
program $P$ hasn’t been run on any of the computer installations from clique),
but corresponding row in $PC$ matrix for $P$ program is non-empty then one
need to use Ridge Regression for the prediction. See Step 4. If the row for
$P$ program is empty then method described in Section 7 is used.

The error of prediction for the algorithm presented above can be estimated by
the equation (1) from Section 2.

###  V-B Matrix decomposition and Tensor decomposition

The problem of programs execution time prediction is very similar to the
problem solved in the recommendation systems. In these systems, there are
usually two types of objects and a relationship between them that can be
quantified. For example, objects can be users and movies, users and books,
users and goods, and so on. The relationship between them is often a rating
that the user puts on a movies, books, or goods. You can build a rating matrix
where the rows (or columns) correspond to movies, books, or goods, and the
columns (or rows) correspond to users. This matrix is often sparse, since
there are a lot of users and objects, and users can’t physically evaluate all
the objects. In the field of recommendation systems, there is a problem is to
determine the ratings of all users for all objects, in other words, you need
to fill in the empty entries in the rating matrix.

Consider using the following analogy: users are computer installations, films,
books or products are programs, and user ratings are execution times. Thus,
computer installations “evaluate” programs and the smaller the rating
(execution time), the better.

The problem of filling empty entries in recommendation systems is solved by
matrix decomposition method [32]. We propose to use the matrix decomposition
technique to solve the problem of programs execution time prediction.

The result of application matrix decomposition techniques to some matrix is
two or more matrices, the product of which gives a matrix that approximates
the original one. For empty entries of the original matrix, i.e. for unknown
values, the product of the matrices gives values that estimate the unknown
values.

$PC$ matrix decomposition allows one to get a vector representations of
programs and computer installations, which have a remarkable property: the
scalar product of the vector representation of the program and the vector
representation of the computer installation is the program execution time on
the computer installation. Thus, the program execution time can be divided
into two components: one is related to the program itself, the other to the
computer installation.

The vector representations of programs and computer installations is
embeddings [39] of programs and computer installations. Embedding techniques
and methods of applying embeddings are very well known in such areas as NLP
[40], topic modeling [41], recommendation systems [42].

The example of matrix decomposition is shown in Fig. 1. In the Rating Matrix
rows correspond to Users, and columns correspond to objects of ratings –
Items. Each entry in the matrix is either empty or the user’s rating. $K$ is a
matrix decomposition parameter that can be changed to achieve the required
accuracy.

The result of matrix decomposition of Rating Matrix of dimension $N$ by $M$ is
two matrices: User Matrix of dimension $N$ by $K$ and Item Matrix of dimension
$K$ by $M$. The User Matrix contains a vector representation of users in rows,
and the Item Matrix contains a vector representation of items in columns. The
product of the user’s vector representation and the item’s vector
representation is the user’s rating for this item. The $K$ parameter affects
not only the accuracy, but also the number of components in the vector
representation.

Figure 1: Example of matrix decomposition.

The following matrix decomposition algorithms were used to solve the problem
of programs execution time prediction: Singular Value Decomposition (SVD)
[43], Alternating Least Square (ALS) [45]. It was also used the tensor
decomposition algorithm [46], based on adaptation of SVD decomposition for
tensors . For tensor decomposition, a $3$-dimensional tensor should be formed
from the original $PC$ matrix. First, the time axis was divided into time
slots of a certain size. Next, the tensor is formed from several $PC^{\prime}$
matrices (the rows – programs, the columns – computer installations), each of
which is associated with a certain time slot. The entries of the $PC^{\prime}$
matrices contain $1$ if the program execution time on the computer
installation hits the corresponding time slot, and $0$ otherwise.

As a result of experimental evaluation, the ALS algorithm was selected as one
with the minimum error (see equation (1) from Section 2.).

###  V-C Ensemble of algorithms

In Sections 5.1 and 5.2, algorithms for the program execution time prediction
on a certain set of computer installations were proposed. These algorithms can
be combined into an ensemble of algorithms to improve the accuracy of the
prediction [47]. In the article averaging ensemble [47] is used, where the
estimation for execution time calculated for each program, as the average
value of the execution time estimations of the following algorithms:

  1. 1.

Ridge Regression – Ridge;

  2. 2.

Computer installations grouping based on the Pearson correlation – Cliques;

  3. 3.

Matrix decomposition by the ALS algorithm – Decomposition.

##  VI Experiments

In this section the results of experimental studies of the algorithms proposed
in Section 5 are presented.

The purposes of the experiments are analysis of the quality of prediction
results:

  1. 1.

based on grouping computer installations by the Pearson correlation;

  2. 2.

based on ALS matrix decomposition algorithm. Selecting the parameter $K$ – the
number of components in the vector representation of programs and computer
installations (see Section 5.2);

  3. 3.

by the ensemble of algorithms;

  4. 4.

by the proposed algorithms with outliers in the source data.

Data

The data sets described in Section 4 were used in experiments. These data sets
were used to form the $PC$ matrix.

The experimental technique

Implemented algorithms [36] were run with different data sets during
experimentation.

###  VI-A Analysis of the quality of prediction based on grouping computer
installations by Pearson correlation

Data set MPIM2007 (see Section 4) with $13$ programs and $396$ computer
installations were used for the experiments.

The algorithm based on grouping computer installations by Pearson correlation
is very sensitive to the presence of outliers in the data, as well as to what
extend the $PC$ matrix is low-density. In order to analyze the quality of
prediction by this algorithm for each entry of the $PC$ matrix, an execution
time prediction was made using information from the remaining entries of this
matrix. Thus, only one value was removed from the matrix, and then the
algorithm described in Section 5.1 was applied. Results are presented in the
Table III.

TABLE III: Results of the prediction based on grouping computer installations by Pearson correlation |  |  | Groups with | Prediction | Average
---|---|---|---|---|---
Exp. | Corr. | Groups | size 1 |  | Error
1 | - | - | - | Regression | 0.25
2 | 0.97 | 46 | 27 | In groups | 0.068
3 | 0.97 | 46 | 27 | In groups + | 0.115
|  |  |  | Regression |

3 experiments (Exp.) were made by following steps:

  1. 1.

Correlation (Corr.) is the threshold for Pearson correlation coefficient
module. It’s necessary to build a graph of computer installations, see Section
5.1;

  2. 2.

Number of groups (Groups) is the number of cliques in graph of computer
installations;

  3. 3.

Number of groups with size 1 (Groups with size 1) is the number of cliques
only with one computer installation;

  4. 4.

Prediction – “Regression” – prediction a value in removed entry by applying of
Ridge Regression to remaining elements in row; “In groups” – prediction in
groups according to the algorithm described in Section 5.1, groups with only
one computer installation aren’t considered in this case; “In groups +
Regression” – prediction execution time in groups with size more than $2$, to
other computer installations used Ridge Regression;

  5. 5.

Average error – after each removing a value from the matrix entry, a
prediction is made, and the prediction error is calculated by equation (1)
from Section 2. After that for all entries the arithmetic mean of prediction
errors for all entries of the matrix is calculated.

According to the results presented in Table III, the algorithm for programs
execution time prediction based on grouping of computer installation by
Pearson correlation gives a prediction error on dense matrices of 11.5% (only
single empty entry in PC matrix). Also the results in Table III shows that the
number of groups with size $1$ is quite large, therefore the prediction error
with and without them differs almost twice: $11.5\%$ and $6.8\%$ respectively.

###  VI-B Analysis of the quality of prediction based on ALS matrix
decomposition algorithm

Figure 2: Ridge regression and ALS matrix decomposition with K=1,2,3,4.

Data set MPIM2007 (see Section 4) with $13$ programs and $396$ computer
installations used for the experiments.

To study the quality of prediction based on ALS matrix decomposition algorithm
$4$ experiments were made: ALS matrix decomposition of $PC$ matrix with
$K=1,2,3,4$. The results of the decomposition were compared with each other,
as well as with Ridge regression algorithm that was chosen as the basic
prediction algorithm.

The Fig. 2 contains graphs of the relationship between the prediction error
and the percentage of empty entry in the $PC$ matrix for Ridge regression and
ALS algorithm with different $K=1,2,3,4$. Prediction error for the algorithms
was calculated using the equation (1) from Section 2. In Fig. 2 $X$-axis is
the percentage of empty entries in the $PC$ matrix (which randomly was removed
from it), $Y$-axis is the prediction error. In the following figures the
axises has the same meaning. Some values from $X$-axis were removed for
readability.

According to the graphs in the Fig. 2 and results from Table III one can make
a conclusion that the Ridge and Cliques algorithms work well on dense matrices
(up to $1\%$ percentage of empty entries). According to the Fig. 2 using the
matrix decomposition technique (by comparison with Ridge and Cliques) with
$K=1$ gives the best results for sparse matrices in which the number of empty
entries is more than $15\%$. Graph for matrix decomposition with $K=1$ is the
lowest in the Fig. 2. Even if $80\%$ of the entries are removed from the $PC$
matrix, it can be filled in so that the execution time will deviate by no more
than $60\%$ from its true value.

Thus, as a result of experiments, we can conclude that the technique of matrix
decomposition with the parameter $K=1$ gives the best solutions when the
percentage of empty entries in the $PC$ matrix is more than $15\%$.

An important advantage of usage the matrix decomposition technique is the
following. The result of applying this technique is vector representations
(embeddings) of programs and computer installations. These embeddings have
dimension $1$! This means that you can enter total ordering on the performance
on a set of computer installations. The less embedding a particular computer
installation has, the less time the program will run on it. This follows from
the fact that for $K=1$, the program execution time on a specific computer
installation is calculated as the product of two numbers: the program
embedding and the embedding of the computer installation.

###  VI-C Analysis of the prediction quality of an ensemble of algorithms

The ensemble of algorithms is described in Section 5.3. The ensemble includes
the following algorithms: Ridge, Cliques and ALS. The simple arithmetic mean
of the prediction results of the Ridge, Cliques and ALS algorithms was
considered. All three data sets – MPIL2007, MPIM2007 and ACCEL_OMP (see
Section 4) were used for the experiments.

The figures contain graphs of the dependence of the prediction error on the
percentage of empty entries in PC matrix for Ridge regression, Cliques, ALS
with K=1, and an ensemble of algorithms.

Figure 3: Results of Ridge, Cliques, Decomposition and an ensemble of
algorithms on MPIL2007 ($1\%$-$24\%$).

The Fig. 3 shows the results for the percentage of empty entries from $1\%$ to
$25\%$. The Fig. 3 shows that the ensemble and Cliques gives the best result
when the percentage of empty entries is up to $4\%$. When this percentage is
increased, the ALS algorithm gives the best result.

Figure 4: Results of Ridge, Cliques, Decomposition and an ensemble of
algorithms on MPIL2007 ($1\%$-$94\%$).

The Fig. 3 shows the results for the percentage of empty entries from $1\%$ to
$94\%$. The Fig. 3 shows that the ALS algorithm gives the best results in all
cases when percentage of empty entries in $PC$ matrix is more than $14\%$.

Figure 5: Results of Ridge, Cliques, Decomposition and an ensemble of
algorithms on MPIM2007 ($1\%$-$25\%$).

In Fig. 5 $X$-axis contains percentage of empty entries in PC matrix from
$1\%$ to $24\%$. The Fig. 5 shows that Ridge regression and Cliques algorithms
give the best results up to $14\%$, and the ensemble of algorithms gives the
best results after $14\%$.

Figure 6: Results of Ridge, Cliques, Decomposition and an ensemble of
algorithms on MPIM2007 ($1\%$-$94\%$).

In Fig. 6 $X$-axis contains percentage of empty entries in PC matrix from
$1\%$ to $94\%$. The Fig. 6 shows that ensemble of algorithms gives the best
results after $14\%$ up to $94\%$.

To test the ensemble of algorithms another experiments was conducted on the
data set ACCEL_OMP with $15$ programs and $25$ computer installations.

Figure 7: Results of Ridge, Cliques, Decomposition and an ensemble of
algorithms on ACCEL_OMP ($1\%$-$92\%$).

Due to small size of $PC$ matrix – $15$ programs, $25$ computer installations
– different algorithms give better results for different percentage of empty
entries in PC matrix. It’s difficult to distinguish the dominance of any one
algorithm. The Fig. 7 shows that first the best result is given by ALS (up to
$19\%$), then ensemble (from $19\%$ to $56\%$), then ALS (from $56\%$ to
$83\%$), then the ensemble again (from $83\%$ to $92\%$).

###  VI-D Analysis of the prediction quality of the proposed algorithms in the
presence of outliers in the source data

The algorithm for program execution time prediction based on grouping computer
installations based on Pearson correlation, like any kind of correlation, is
sensitive to outliers in the source data. Therefore, two more experiments were
conducted to test the stability of each algorithm to outliers.

For the experiments, the PC matrix was formed, after that $10\%$ of its
entries were randomly selected and multiplied:

  * •

In the first experiment on a random number in the interval $(0,4)$;

  * •

In the second experiment on a random number in the interval $(0,10)$.

Figure 8: Results of Ridge, Cliques, Decomposition and an ensemble of
algorithms on MPIM2007 ($1\%$-$92\%$), $10\%$ of outliers, coefficient from
$0$ to $4$. Figure 9: Results of Ridge, Cliques, Decomposition and an ensemble
of algorithms on MPIM2007 ($1\%$-$92\%$), $10\%$ of outliers, coefficient from
$0$ to $10$.

The Fig. 8 and Fig. 9 shows the results of a prediction based on data with
outliers. $X$-axis contains percentage of empty entries in PC matrix. $Y$-axis
contains prediction errors. As you can see from the Fig. 8 and Fig. 9 the
ensemble of algorithms is almost always better than all other algorithms. The
Cliques algorithm also gives a good result up to $19\%$ of the percentage of
empty entries.

###  VI-E Conclusions from the experiments

Experiments in the Sections 6.1-6.4 showed that an ensemble of algorithms with
a small percentage of empty entries in “Program-Computer” matrix up to 4% and
middle percentage of empty entries – from 19% to 56% – makes a better
prediction compared to all other algorithms. Also, as experiment in Section
6.4 showed, the ensemble and ALS show good results even in the presence of
outliers in the source data sets. When the percentage of empty entries is much
more than 50%, the best prediction is often made by the ALS algorithm.

Thus, for dense matrices, it is better to use an ensemble of algorithms, for
sparse ones – matrix decomposition, in particular the ALS algorithm.

##  VII Execution time prediction application

Based on the research presented above the following algorithm was constructed
for solving the problem of efficient resource allocation in a cloud computing
environment over a heterogeneous network of computers:

  1. 1.

Create a $PC$ matrix based on the history of finished program executions;

  2. 2.

Fill in empty entries in $PC$ matrix using one of the algorithms described in
the article (see Section 5);

  3. 3.

Knowing the estimated programs execution times on computer installations apply
one of the algorithms below:

    1. (a)

Greedy algorithm. Send the current program to the computer installation where
it has the minimum execution time. If this program hasn’t been executed
anywhere else, it can be placed on the most powerful computer installation.
Performance can be determined by the embedding of the computer installation
(see Section 5);

    2. (b)

Build a schedule. If one plan the calculation for multiple programs, one can
create a schedule for running Programs that minimizes the total execution time
[12].

The experiments in Section 6 used “fixed” data. A number of programs and a
number of computer installations were fixed. But one can not fix them, and
change the $PC$ matrix in the process (constantly). So there are two ways to
get data for the prediction algorithm:

  1. 1.

Take an arbitrary set of benchmarks from different subject areas and different
“behavior” with the processor, memory, and network. Start them on all computer
installations from the specified set. Build a $PC$ matrix, make a matrix
decomposition, and get vector representations of programs and computer
installations. Vector representations of computer installations can be used to
rank computer installations by inverse “performance” if the vector
representation has dimension $1$. For new programs that don’t have execution
information in the $PC$ matrix, you can predict execution time using only on
the ranks of the computer installations;

  2. 2.

When choosing a computer installation to start a program, you need to evaluate
the program execution time, and not just rely on the ranks of the computer
installations. This is important because the behavior of the same program can
be completely unpredictable on different computer installations. This is due
to the fact that different computer installations can perform certain
functions, operations, and commands that are present in the program in
different ways. Therefore, to take into account the “behavior” of the program,
you need to constantly change the matrix: add new rows for new programs, new
rows for old programs, but with new arguments, new columns for new computer
installations.

Set of benchmarks from the website spec.org [48] can be used to build a $PC$
matrix.

In the experiments presented in Section 6, each program had only one set of
arguments and only one row in the $PC$ matrix corresponded to it. Such data
was considered to for simplicity of studying the quality of the created
algorithms. However, the developed algorithms can be easily applied for the
case when each program was run several times. There are at least $2$ ways to
apply the described algorithms in this case:

  1. 1.

Form a $PC$ matrix as described in Section 2. In this case, a single program
will correspond to a group of strings. Then the algorithms proposed in this
article in Section 5 are applied;

  2. 2.

In case of several program executions, you can first make predictions based on
the history – predictions of the program execution time within a single
computer installation [22], [23], [24], [25]. This way, you can fill in the
empty entries in the columns. Then the algorithms proposed in this article in
Section 5 are applied to fill the remaining empty entries.

##  VIII Conclusion

This article describes algorithms for predicting the Programs execution time
on the federation’s computer installations that meet the requirements of the
virtual infrastructure. Existing algorithms for predicting program execution
time use the history of executions to make a prediction. The main drawback of
existing algorithms is that the algorithms work only within a single computer
installation, i.e. to predict the execution time on a certain computer
installation, you need more than one start of the program on this computer
installation – that is, you can’t do without the history.

In this paper, a new approach that allows to predict the Program execution
time on a certain computer installation was proposed, even when it was not
executed on it. Two algorithms were constructed and analyzed: an algorithm
based on computer installations grouping based on the Pearson correlation
coefficient (for cases when the Program was executed on almost all federates
in the federation) and an algorithm based on matrix decomposition technique,
which allows to obtain vector representations (embeddings) of the Program and
computer installations (for cases when the Program was executed only on a
small subset of federates). Matrix decomposition, which was used to predict
the Programs execution times, is actively used in the field of recommendation
systems, for example, SVD decomposition [43], ALS [45] and tensor
decomposition [46]. In the article, for the first time, it is proposed to use
this technique for the problem of predicting the programs execution time on an
set of computer installations.

In addition, an ensemble of algorithms consisting of Ridge regression,
grouping computer installations based on Pearson correlation coefficient, and
matrix decomposition was proposed. It is shown that the ensemble of algorithms
is more resistant to outliers than other algorithms and gives the best results
on dense matrices.

Embeddings of programs and computer installations were obtained as a result of
matrix decomposition of the matrix ”Programs-Computers”, where value in each
entry is the program execution time on the computer installation. This
approach to obtaining embeddings is widely used in many areas, for example, in
recommendation systems [42], in thematic modeling [41]. This approach was
successfully applied in solving the problem. The approach proved to be
universal, as similar results were obtained for MPI benchmarks and OpenMP
benchmarks (see experiments in Section 6).

It is important to emphasize that the proposed approach to predicting program
execution time requires a minimal set of data about the program, which is
usually available on all modern computer installations. Another important
advantage of the proposed approach is that a result of matrix decomposition is
the embeddings of Programs and computer installations of dimension 1. This
fact allows one set up the total order as on a set computer installations as
on a set of programs what significantly help to properly select the computer
installation with effective execution time.

The final remark is as a hypothesis we state that execution time prediction
techniques proposed in this article one can apply not only to MPI programs.

## References

  * [1] Hwang, T. (2017, March). NSF GENI cloud enabled architecture for distributed scientific computing. In 2017 IEEE Aerospace Conference (pp. 1-8). IEEE
  * [2] Baldin, Ilya and Nikolich, Anita and Griffioen, James and Monga, Indermohan and Wang, Kuang-Ching and Lehman, Tom and Ruth, Paul. ”FABRIC: A National-ScaleProgrammable ExperimentalNetwork Infrastructure,” IEEE Internet Computing, v.23, 2020
  * [3] V. Antonenko, R. Smeliansky, I. Baldin, Y. Izhvanov and Y. Gugel, ”Towards SDI-bases Infrastructure for supporting science in Russia,” 2014 International Science and Technology Conference (Modern Networking Technologies) (MoNeTeC), Moscow, 2014, pp. 1-7, doi: 10.1109/MoNeTeC.2014.6995576.
  * [4] Dewar, R. G., MacKinnon, L. M., Pooley, R. J., Smith, A. D., Smith, M. J., & Wilcox, P. A. (2002). The OPHELIA Project: Supporting Software Development in a Distributed Environment. In ICWI (pp. 568-571).
  * [5] European Commission, “Federation for FIRE, projects from FP7 programme” [Online]. Available: https://cordis.europa.eu/project/rcn/ 105823/ [Accessed: 29-Apr-2019]
  * [6] European Commission, “Federation for FIRE plus, projects from H2020 programme” [Online]. Available: https://cordis.europa.eu/project/rcn/ 207030/ [Accessed: 29-Apr-2019]
  * [7] Bal, H. E., Bhoedjang, R., Hofman, R., Jacobs, C., Langendoen, K., Rühl, T., & Kaashoek, M. F. (1998). Performance evaluation of the Orca shared-object system. ACM Transactions on Computer Systems (TOCS), 16(1), 1-40.
  * [8] Brun, R., Urban, L., Carminati, F., Giani, S., Maire, M., McPherson, A., Patrick, G. (1993). GEANT: detector description and simulation tool (No. CERN-W-5013). CERN
  * [9] Antonenko, V., Petrov, I., Smeliansky, R.L., Huang, Z., Chen, M., Cao, D., & Chen, X. (2019). MC2E: META-CLOUD COMPUTING ENVIRONMENT FOR HPC
  * [10] Ridge Regression[Online] Available: https://ncss-wpengine.netdna-ssl.com/wp-content/themes/ncss/pdf/Procedures/NCSS/Ridge_Regression.pdf [Accessed: 21-Jun-2020]
  * [11] Random forest algorithm [Online] Available: https://www.stat.berkeley.edu/ breiman/RandomForests/cc_home.htm [Accessed: 21-Jun-2020]
  * [12] MR Gary and DS Johnson, Computers and Intractability: A Guide to the Theory of NP-CompletenessJanuary 1979
  * [13] González, Marcos Tulio Amarís. Performance prediction of application executed on GPUs using a simple analytical model and machine learning techniques (2018) [Online] Available: http://doi.org/10.11606/t.45.2018.tde-06092018-213258 [Accessed: 21-Jun-2020]
  * [14] A. Snavely, L. Carrington, N. Wolter, J. Labarta, R. Badia and A. Purkayastha, ”A Framework for Performance Modeling and Prediction,” SC ’02: Proceedings of the 2002 ACM/IEEE Conference on Supercomputing, Baltimore, MD, USA, 2002, pp. 21-21, doi: 10.1109/SC.2002.10004
  * [15] Nadeem, F., Alghazzawi, D., Mashat, A., Fakeeh, K., Almalaise, A., & Hagras, H. (2017). Modeling and predicting execution time of scientific workflows in the Grid using radial basis function neural network. Cluster Computing, 20(3), 2805–2819. doi:10.1007/s10586-017-1018-x
  * [16] Brehm, J., Worley, P.H.: Performance prediction for complex parallel applications. In: Proceedings of the 11th International Parallel Processing Symposium, 1997, pp. 187–191. IEEE (1997)
  * [17] Bacigalupo, D.A., et al.: An investigation into the application of different performance prediction methods to distributed enterprise applications. J. Supercomput. 34(2), 93–111 (2005)
  * [18] Yero, E.J.H., Henriques, M.A.A.: Contention-sensitive static performance prediction for parallel distributed applications. Perform. Eval. 63(4), 265–277 (2006)
  * [19] Barnes, B.J., Rountree, B., Lowenthal, D.K., Reeves, J., de Supinski, B., Schulz, M.: A regression-based approach to scalability prediction. In: Proceedings of the 22nd Annual International Conference on Supercomputing, ICS ’08, pp. 368–377. ACM (2008)
  * [20] Wu, Q., Datla, V.V.: On performance modeling and prediction in support of scientific workflow optimization. In: Proceedings of the 2011 IEEE World Congress on Services, SERVICES ’11, pp. 161– 168. IEEE Computer Society (2011)
  * [21] Wu, X., Taylor, V., Paris, J.: A web-based prophesy automated performance modeling system. In: Proceedings of the International Conference onWeb Technologies, Applications and Services (WTAS2006), pp. 17–19 (2006)
  * [22] Gibbons, R.: A historical application profiler for use by parallel schedulers. In: Proceedings of the Job Scheduling Strategies for Parallel Processing. Springer (1997)
  * [23] Kapadia, N.H., Fortes, J.A., Brodley, C.E.: Predictive application performance modeling in a computational grid environment. In: Proceedings of the Eighth International Symposium on High Performance Distributed Computing, pp. 47–54. IEEE (1999)
  * [24] Li, H., Groep, D., Templon, J., Wolters, L.: Predicting job start times on clusters. In: CCGRID ’04: Proceedings of the 2004 IEEE International Symposium on Cluster Computing and the Grid, pp. 301–308 (2004)
  * [25] Mohr, B., Wolf, F.: Kojak—a tool set for automatic performance analysis of parallel programs. In: Euro-Par 2003 Parallel Processing, pp. 1301–1304. Springer (2003)
  * [26] Iverson, M.A., Özgüner, F., Potter, L.: Statistical prediction of task execution times through analytic benchmarking for scheduling in a heterogeneous environment. IEEE Trans. Comput. 48(12), 1374–1379 (1999)
  * [27] Liu, X., Chen, J., Liu, K., Yang, Y.: Forecasting duration intervals of scientific workflow activities based on time-series patterns. In: Proceedings of the IEEE Fourth International Conference on eScience, 2008, eScience ’08, pp. 23–30 (2008)
  * [28] Ipek, E., de Supinski, B.R., Schulz, M., McKee, S.A.: An approach to performance prediction for parallel applications. In: Proceedings of the 11th International Euro-Par Conference on Parallel Processing, Euro-Par’05, pp. 196–205. Springer, Berlin, Heidelberg (2005)
  * [29] Lee, B.C., Brooks, D.M., de Supinski, B.R., Schulz, M., Singh, K., McKee, S.A.: Methods of inference and learning for performance modeling of parallel applications. In: Proceedings of the 12th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming, PPoPP ’07, pp. 249–258. ACM (2007)
  * [30] Yang, L.T., Ma, X., Mueller, F.: Cross-platform performance prediction of parallel applications using partial execution. In: Proceedings of the 2005 ACM/IEEE Conference on Supercomputing, SC ’05, pp. 40–44. IEEE Computer Society (2005)
  * [31] Pearson’s Correlation Coefficient [Online] Available: https://link.springer.com/referenceworkentry/10.1007
  * [32] Cheng CM., Jin XQ. (2018) Matrix Decomposition. In: Alhajj R., Rokne J. (eds) Encyclopedia of Social Network Analysis and Mining. Springer, New York, NY
  * [33] MPI2007 datasets [Online] Available: https://spec.org/mpi2007/results/mpi2007.html [Accessed: 24-Jun-2020]
  * [34] Accel OpenMP dataset [Online] Available: https://spec.org/accel/results/accel.html [Accessed: 24-Jun-2020]
  * [35] Update history for MPI2007 [Online] Available: https://spec.org/mpi2007/results/ [Accessed: 24-Jun-2020]
  * [36] MPI program execution time prediction algorithms from the article. https://github.com/andxeg/CloudCom_2020_Embeddings
  * [37] Skiena, S. S. ”Clique and Independent Set” and ”Clique.” §6.2.3 and 8.5.1 in The Algorithm Design Manual. New York: Springer-Verlag, pp. 144 and 312-314, 1997
  * [38] Johnston, H.C. Cliques of a graph-variations on the Bron-Kerbosch algorithm. International Journal of Computer and Information Sciences 5, 209–238 (1976). https://doi.org/10.1007/BF00991836
  * [39] Add link to embedding definition
  * [40] Li Y., Yang T. (2018) Word Embedding for Understanding Natural Language: A Survey. In: Srinivasan S. (eds) Guide to Big Data Applications. Studies in Big Data, vol 26. Springer, Cham
  * [41] Aggarwal C.C. (2018) Matrix Factorization and Topic Modeling. In: Machine Learning for Text. Springer, Cham
  * [42] Recommender system. [Online] Available: https://link.springer.com/referenceworkentry/10.1007
  * [43] Helmke U., Moore J.B. (1994) Singular Value Decomposition. In: Optimization and Dynamical Systems. Communications and Control Engineering. Springer, London
  * [44] Pierre Comon, Xavier Luciani, and Andre LF De Almeida (2009) “Tensor decompositions, alternating least squares and other tales”. Journal of Chemometrics 23.7-8, pp.393-205
  * [45] Gábor Takács and Domonkos Tikk. 2012. Alternating least squares for personalized ranking. In Proceedings of the sixth ACM conference on Recommender systems (RecSys ’12). Association for Computing Machinery, New York, NY, USA, 83–90. DOI:https://doi.org/10.1145/2365952.2365972
  * [46] Evgeny Frolov and Ivan Oseledets. 2016. Fifty Shades of Ratings: How to Benefit from a Negative Feedback in Top-N Recommendations Tasks. In Proceedings of the 10th ACM Conference on Recommender Systems (RecSys ’16). Association for Computing Machinery, New York, NY, USA, 91–98. DOI:https://doi.org/10.1145/2959100.2959170
  * [47] Džeroski S., Panov P., Ženko B. (2009) Machine Learning, Ensemble Methods in. In: Meyers R. (eds) Encyclopedia of Complexity and Systems Science. Springer, New York, NY .RIS
  * [48] MPI benchmarks [Online] Available: https://spec.org/mpi2007/benchmarks/ [Accessed: 21-Jun-2020]