# ParaGraph: Weighted Graph Representation for Performance Optimization of HPC
Kernels

Ali TehraniJamsaz1, Alok Mishra2, Akash Dutta1, Abid M. Malik3, Barbara
Chapman2, Ali Jannesari1
1Iowa State University, Ames, Iowa, USA {tehrani, adutta,
<EMAIL_ADDRESS>2Stony Brook University, Stony Brook, New York, USA
{almishra<EMAIL_ADDRESS>3Brookhaven National Laboratory,
Upton, New York, USA<EMAIL_ADDRESS>
###### Abstract

GPU-based HPC clusters are attracting more scientific application developers
due to their extensive parallelism and energy efficiency. In order to achieve
portability among a variety of multi/many core architectures, a popular choice
for an application developer is to utilize directive-based parallel
programming models, such as OpenMP. However, even with OpenMP, the developer
must choose from among many strategies for exploiting a GPU or a CPU.
Recently, Machine Learning (ML) approaches have brought significant advances
in the optimizations of HPC applications. To this end, several ways have been
proposed to represent applications’ characteristics for ML models. However,
the available techniques fail to capture features that are crucial for
exposing parallelism. In this paper, we introduce a new graph-based program
representation for parallel applications that extends the Abstract Syntax Tree
to represent control and data flow information. The originality of this work
lies in the addition of new edges exploiting the implicit ordering and parent-
child relationships in ASTs, as well as the introduction of edge weights to
account for loop and condition information. We evaluate our proposed
representation by training a Graph Neural Network (GNN) to predict the runtime
of an OpenMP code region across CPUs and GPUs. Various transformations
utilizing collapse and data transfer between the CPU and GPU are used to
construct the dataset. The predicted runtime of the model is used to determine
which transformation provides the best performance. Results show that our
approach is indeed effective and has normalized RMSE as low as $4\times
10^{-3}$ to at most $1\times 10^{-2}$ in its runtime predictions.

###### Index Terms:

OpenMP, HPC, offloading, program representation

##  I Introduction

Over the years, the increase in the number of on-chip cores has led to
significant improvement in the performance of parallel code. To exploit this
increased computation capacity, applications have been readily modified using
Pthreads or OpenMP constructs. In the last decade, General Purpose Graphic
Processing Units (GPGPU) have gained popularity. They can handle massive data
parallelism with low power consumption. As a result, most HPC platforms are
now coupled with GPGPU(s). HPC platforms will continue to support more
accelerators, but given the difficulties in utilizing and configuring even one
accelerator, using and configuring multiple heterogeneous accelerators will
become increasingly difficult in the future. On the other hand, utilizing GPUs
effectively imposes challenges that require re-engineering the code and
applications. The full computing power of GPUs may not be utilized if they are
not used properly. In addition, improper data transfer between the host (CPU)
and device (GPU) can result in significant overheads. It is exceptionally
challenging and burdensome for developers to create applications for extremely
heterogeneous platforms with multiple devices. The recent emergence of tools
and programming models aims to automate this process of application adaptation
to heterogeneous platforms. One of the most popular parallel programming
models, OpenMP [1], aims to make the process of developing parallel programs
that can run on different architectures simpler. Despite this, optimizing a
code to use the OpenMP directives correctly is still a tedious task for large
and complex applications.

Recent advances in Deep Learning (DL) have enabled researchers to apply DL to
a wide range of software engineering jobs and challenges, including code
comment generation, compiler optimizations, and heterogeneous platform
adaptation. Since source code cannot be directly fed into DL models, we need a
suitable representation for applications to serve as input to various DL
models. In this paper, we propose $ParaGraph$, a graph-based program
representation that aims to expose critical characteristics (e.g., the number
of iterations in a loop) of HPC applications. To the best of our knowledge,
existing program representations are not designed to expose and represent the
characteristics of parallel HPC applications. These program representations
usually target serial code and do not incorporate characteristics of a
parallel code. As a result, DL models built on top of these representations
cannot model features inherent to parallel codes.

Previous works, such as COMPOFF [2], have relied on feature engineering to
overcome these shortcomings. However, feature engineering requires expert
knowledge. For a fast-evolving field such as HPC, always relying on expert
intervention for such feature engineering is not realistic. There is a need
for an adaptable approach that can automatically extract such features. Our
proposed approach aims to address this and uses Graph Neural Networks (GNNs)
to model the code graphs generated by $ParaGraph$. Experimental results show
that our model can predict the runtime of HPC kernels with a very low error
rate (at most $1\times 10^{-2}$ in terms of normalized RMSE), confirming the
efficacy of our strategic approach.

In the rest of the paper, we first discuss some background work in Section II
on program representation, OpenMP GPU offloading, DL in compilers, OpenMP
Advisor (a tool we used to generate various kernels for our work), and other
related work. Then in Section III, we define $ParaGraph$. Section IV covers
the experiments carried out in this paper, and Section V discusses the result.
Finally, we conclude our work with discussions of future plans in Section VI.

##  II Background and Other Related Work

In this section, related works and some background for program representation,
OpenMP parallelism, and its tools are discussed.

###  II-A Program Representation

Recently, with breakthroughs in Machine Learning (ML) and specifically Deep
Learning (DL), researchers have been applying data-driven approaches to
various software engineering tasks, and challenges, ranging from code comment
generation [3] to compiler optimizations [4]. However, the source code can not
simply be fed to DL models as DL models do not accept strings as input.
Therefore, applications need to be presented in a way usable by the models.
Some works have presented programs as a sequence of tokens [5, 6]. While this
kind of representation has worked well for some downstream tasks, presenting
programs as a sequence of tokens discards critical syntactical and semantic
characteristics of programs which are essential for deep learning models to
reason.

Recently, to better present programs’ semantics and syntactic information,
graph-based representations have been proposed. Allamanis et al. [7] proposed
enhanced ASTs and added new edges to present control flow information. They
applied the augmented AST to two tasks: variable misuse detection and variable
name suggestion. Inspired by their work, we also augment the AST by adding new
edges. However, we incorporate additional edges to better represent $loops$
and $if$ conditions.

Cummins et al. in [4] proposed a lower-level graph representation based on
LLVM intermediate representation for solving compiler optimization tasks.These
representations are effective for downstream tasks, such as simple algorithm
classification. To the best of our knowledge, there does not exist a
representation tailored toward representing the characteristics of parallel
loops and if-statements.

###  II-B OpenMP GPU Offloading

When different architecture supports either a different native language or
various optimization techniques for the same language, program portability
becomes a serious problem. Regrettably, programming languages and compilers
are still far from being able to handle portability on their own. Program
portability should be the primary priority for developers who have to deal
with multiple device memory accesses on the same data or someone whose
programs must work effectively on systems with various node architectures,
such as many-core vs. GPU-accelerated nodes.

Utilizing a directive-based programming paradigm, such as OpenMP, the de-facto
standard for parallel programming in C/C++ and Fortran, is one approach to
ensuring portability across several architectures. OpenMP intends to move to
extremely heterogeneous architectures [8] and has supported GPU offloading
from specification 4.0. Unfortunately, even with OpenMP, optimizing large-
scale applications remains a challenging task.

parallel forparallel forparallel forparallel fordistributetarget teamstarget
region ends

Figure 1: OpenMP target teams distribute parallel for.

Since GPUs are very parallel machines, programmers should take full advantage
of this fact to get the most performance possible. The “omp parallel for”
pragma alone will only parallelize a code for CPUs; it won’t offload the
computation to a GPU. We expect a high level of coarse grain parallelism on a
platform like GPU. The amount of parallelism that a GPU can use is constrained
by this design. Figure 1 illustrates how OpenMP teams and distribute
directives create additional levels of parallelism. At the start of a target
region, only one team and one member thread are active. The teams distribute
directive distributes the full loop iteration space among all available teams.
Additionally, in the case of nested parallel loops, we can utilize the
parallel for directive on them to further distribute the iterations of the
nested loops among threads within a team. We utilize the combined directive
“teams distribute parallel for” to distribute the iteration space of one loop
among teams and threads inside a team when there is just one level of parallel
loops or when the outer loop has sufficient parallelism. teams are used to
group threads, and distribute enables a team group to be scheduled to run in a
loop. Teams resemble CUDA threadblocks on NVIDIA GPUs, as there is no
synchronization primitive to function as a barrier between threads from
different teams. And the threads within each team are parallelized typically
using parallel for.

There are several frameworks being worked on right now that will automatically
assist application developers in handling severe heterogeneity. These
frameworks necessitate analytical cost models, which will help them select the
kind of optimization the application requires. Hand-tuned cost functions are
extensively employed currently; however, calculating optimization costs
requires a deeper understanding of the underlying hardware. Despite its
efficacy, manually building a cost model for a single architecture can take
months. Since cost functions are crucial and manual tuning is time-consuming,
compiler engineers are investigating Machine, and Deep Learning approaches as
a means of automating this process.

###  II-C ML in Compiler

Recently, a lot of research has been done on how to use learning-based
techniques in compilation as well. Early work exploiting ML in compilers
primarily explored its use to help optimize sequential programs. However, its
application to the task of optimizing parallel programs has attracted
significant attention in the past decade due to the prevalence of multi-core
platforms and, more recently, heterogeneous systems [9, 10]. Mirka et.al. [11]
devised a decision-tree-based approach to predict the scheduling policy for an
OpenMP parallel region. Also, Denoyelle et.al. [12] uses machine learning
techniques to optimize OpenMP programs for scheduling policies and the number
of threads. Alcaraz et.al. [13] uses neural networks and performance counters
to predict the number of threads. Dutta et.al. [14] uses an LLVM-IR based
graph representation and performance counter to train a GNN model to predict
number of thread, chunk size and scheduling policy for OpenMP loops. Tree and
graph-based features have also been used by Malik et.al. [15], who present a
unique graph-based approach for feature representation. Learning-based
techniques were used to build classifiers to determine whether to offload
OpenCL code [16] and to select a clock frequency at which the processor should
operate [17]. A high level of accuracy was reported; however, the benefits
could not be quantified as the work did not attempt to generate modified code.
They also explored regression techniques to build curve fitting models to
search for the sweet spot for work partitioning among processors [18] or a
trade-off of energy and performance [19].

The results of prior efforts from applying ML on compiler optimizations are
encouraging. However, new feature engineering practices need to be explored
that can help learning-based methods to learn more about a code and its
computational needs. In COMPOFF [2], the authors provide a proof of concept
for using an Artificial Neural Network model to make better decisions in
offloading a kernel to a GPU using OpenMP. They utilize a fully-connected
feed-forward network, also referred to as multi-layer perceptrons (MLPs),
which are effectively stacked layers of linear regression.

###  II-D OpenMP Advisor

OpenMP Advisor [20] is a first-of-its-kind compiler tool that enables OpenMP
code offloading to a GPU via Machine Learning. The Advisor is divided into
three major modules: Kernel Analysis, Cost Model, and Code Transformation. The
Kernel Analysis module recommends various variants for a given application,
and the Code Transformation module generates codes for those variants. In our
work, we use the code transformation module of OpenMP Advisor to generate
various kernels for training our model. Through the use of COMPOFF, a machine
learning-based cost model, this tool can identify the kernels that are most
suitable for offloading by predicting their runtime. But COMPOFF has some
limitations. It requires figuring out how many operations are contained within
a kernel, which is a challenging task in and of itself. The Advisor’s current
functionality is only restricted to GPUs, though it has the potential to be
expanded to support additional OpenMP-capable hardware. One technique for
expanding to other devices is to develop a model that works outside GPUs as
well. Extending COMPOFF to work beyond GPUs is a daunting task, and there is a
need for a new cost model that can bridge this gap.

###  II-E Other Related Works

A recent tool called BLISS [21] allows for auto-tuning parallel applications
without the use of instrumentation or domain-specific knowledge. This auto-
tuner develops a variety of lightweight models with the aid of Bayesian
optimization that can compete with the output quality of a complex,
heavyweight Machine Learning (ML) model. In the same context, Bayesian
optimization is frequently used by other autotuners, including BOHB [22],
HiPerBOt [23], GPTune [24] and ytopt [25]. Unfortunately, due to their
expensive evaluation functions, tuning large-scale HPC systems, which are
becoming increasingly heterogeneous, complex, and expensive to evaluate, is
still quite challenging. Nonetheless, most of the autotuners mentioned above
are online, meaning that even in the inference time, they execute applications
iteratively to tune them. Therefore, in terms of complex heterogeneous large-
scale HPC systems, these tuners will be very costly to apply. In contrast,
$ParaGraph$ is an offline approach and does not require the execution of
applications in the inference phase. COMPOFF [2] provides a novel portable
cost model that statically calculates the cost of OpenMP offloading on various
architectures. COMPOFF is also an offline approach. In this work, we compare
$ParaGraph$ to COMPOFF since both can predict the runtime of a kernel on a GPU
statically.

##  III ParaGraph

int x;// ...x =
50;CompoundStmtBinaryOperatorImplicitCastExprIntegerLiteralVarDecl. .
.DeclRefExprCompoundStmtBinaryOperatorImplicitCastExprIntegerLiteralVarDecl. .
.DeclRefExpr1111NextTokenNextSiblingRefif ( x > 50) {/*...*/}else
{/*...*/}CompoundStmtIfStmtBinaryOperatorCompoundStmtCompoundStmtCompoundStmtIfStmtBinaryOperatorCompoundStmtCompoundStmt50502525ConTrueConFalsefor(int
i=0; i<50; i++) {
/*...*/}CompoundStmtForStmtDeclStmtBinaryOperatorCompoundStmtUnaryStmtCompoundStmtForStmtDeclStmtBinaryOperatorCompoundStmtUnaryStmt11505050ForExecForExecForNextForNext
AST Edge NextToken Reference NextSibling ConFalse ConTrue ForExec ForNext

Figure 2: Modification to AST to create Augmented AST for ParaGraph

In this section, we present $ParaGraph$, a weighted graph representation of
programs that captures characteristics related to HPC kernels. Abstract Syntax
Tree serves as a base for $ParaGraph$. In fact, additional edges and
attributes are introduced to AST to construct $ParaGraph$. Some previous works
have proposed various ways of augmenting ASTs to present semantic information
[26, 7]. Although these AST-based program representations are effective in
training deep learning models for downstream tasks such as variable misuse
detection, not enough attention has been given to better represent regions of
programs that have a high impact on the performance of the application, such
as $loops$ and $if$ condition. In $ParaGraph$ we try to address these
shortcomings. As stated, $ParaGraph$ is based on AST; however, it includes new
edges to present characteristics such as control flow or order of children of
a node. Moreover, $loops$ are represented by special edges conveying the order
in which the loop condition and its body execute iteratively. Additionally, we
add weights to the edges to present how many times each node will be utilized
during the execution of a program. It is worth noting that all these
augmentations are applied statically; therefore, they do not cause significant
overhead. $ParaGraph$ can be directly utilized by Graph Neural Network models,
and our experimental results show that, in performance optimization, this
representation is quite effective and outperforms the state-of-the-art
approach.

###  III-A ParaGraph Construction

####  III-A1 Abstract Syntax Tree

As stated, the base of $ParaGraph$ is Abstract Syntax Tree. Typically, any
kind of compiler can be used to generate AST. In this study, we use Clang, a
popular C/C++ front-end, to parse and compile C/C++ programs. ASTs contain two
types of nodes: non-terminal and terminal. Non-terminal nodes are often called
syntax nodes, whereas terminal nodes are called syntax tokens. The relation
between nodes in AST is a parent-child relation. Once an AST is retrieved, we
add additional edges and attributes depicting information related to control
flow and data flow. The following subsections will provide more details on
these augmentations.

####  III-A2 Augmenting Abstract Syntax Tree

ASTs typically provide structural and syntactical information of a program.
Although this information can be useful for neural networks to learn the
characteristics of programs to some extent, it has been shown that adding
additional attributes and information, especially from a compiler point of
view, boosts the learning capabilities of deep learning models. Inspired by
existing works [26, 7], we introduce the following additional edges to AST
(shown in Figure 2) to represent information related to the control and data
flow of the programs.

  * –

NextToken:  By default there is no order imposed among the syntax tokens.
However, from a compiler’s perspective, syntax tokens have an order. The order
among the syntax tokens is from far left token to the far right one. To
present this information in a graph, we introduce a new edge type called
NextToken. NextToken connects each syntax token to the syntax token on its
right side. In other words, it connects each syntax token to the next syntax
token. In this way, an order is imposed on syntax tokens. This edge type is
shown in orange in Figure 2 (the AST on the left).

  * –

NextSib:  As mentioned above, edges in ASTs show a parent-child relationship.
Compilers assume an order among children nodes. For example, a binary operator
such as division has two children. In an AST, the left child is always the
numerator, and the right child is the denominator. Therefore, it is necessary
to show the order of the children. NextSib presents this information. NextSib
connects each syntax node to its sibling on the right (the blue edge in Figure
2.

  * –

Ref:  In clang’s AST, referenced variables are presented by DecRefExpr nodes.
These nodes are terminal and do not have children. During graph construction,
we introduce Ref edges (shown in pink in Figure 2) connecting a DecRefExpr
node to where the corresponding variable is defined. This conveys information
where a declared variable is used throughout the graph.

  * –

ForNext, ForExec:  Loops are typically shown as ForStmt nodes in Clang’s AST.
ForStmt nodes usually have four children. First child is related to
initializing the loop counter, whereas the second child is about checking the
loop’s condition, deciding whether the next iteration should be executed or
the loop has ended. Third child is CompoundStmt which presents the body of a
loop. Lastly, the last child is a modifier that modifies the loop counter such
as incrementing it by one. The relationship between these four children
exposes critical characteristics of loops. These characteristics are known by
compilers but are not presented explicitly in the ASTs. We add ForExec edges,
which connect the first child (counter initialization) to the second child
(loop condition) and also connect the second child to the third child (loop
body) as shown in the right side of Figure 2. In fact, ForExec edges show the
flow of executing the next iteration of the loop. We connect the third child
(Loop’s body) of the ForStmt node to the fourth child (loop counter modifier)
via a new type of edge called ForNext. Basically, ForNext represents the
information related to whether the next iteration of the loop needs to be
executed or not, whereas ForExec represents the next execution of the loop.
After the fourth child (loop counter modifier) comes the second child (loop
condition), which checks if next iteration is going to happen or if the loop
has ended. Therefore, the fourth child is connected to the second child with a
ForNext edge as shown in Figure 2.

  * –

ConTrue, ConFalse:  If statements in the ASTs are shown by IfStmt nodes.
IfStmt nodes typically have three children. First child presents the
comparison condition; the second child is the body of the if condition, and
the third child is the body of the else part. To present this information, we
connect the first child to the second child through a new type of edge called
ConTrue to show that the flow moves to the second child if the if statement is
true. On the other hand, the first child is connected to the third child via a
ConFalse edge. This edge shows that if the condition of if statement is not
met, the flow moves to the third child.

  * –

Child: Child edges are normal AST edges that present a parent-child
relationship between the nodes in AST.

ASTInputAST Variant Generator ASTASTASTASTASTASTVariantsAST Compiler Frontend
ASTAST ParaGraph Generator
ASTAST11NN111ASTASTASTASTASTASTASTParaGraphASTParaGraph GenerationAST Compile
and Build AST Runtime Measurement Module ASTASTTime
(ms)AST746846AST9984695AST916548ASTMeasuring Execution Time of
ApplicationsASTASTASTASTASTASTASTASTASTASTASTASTDatasetAST GNN model
ASTPredictionASTPredictAST T r a i n ASTASTTrainAST

Figure 3: Workflow of ParaGraph

####  III-A3 Weighted Edges

In the previous section, we explained how we augment the AST by introducing
new types of edges to better present programs. However, there is still some
information that is missing in our representation. For instance, whenever we
encounter a loop in a program, the loop’s body could be executed multiple
times. Or the branches of an if statement will not execute the same number of
times. To present this information in our graph structure, we augment the
edges by adding weights. Weights of edges are calculated based on the region
and edge type. Weights are added only to Child edges since these are AST
edges, and a compiler uses the AST nodes and child edges to transform an AST
to lower-level for machine execution. For illustration purpose, some edge
weights are shown in Figure 2. We add the weights as follows:

  * –

Default edge weight: By default, we add a weight of 1 to edges to represent
the fact that each statement in code executes one time only, and once a
statement executes, the control moves to the next statement.

  * –

Loops: Despite the typical statements that execute only once, statements in
the body of a loop will be executed several times depending on the number of
iterations. To expose the number of iterations in our graph representation, we
first observe the number of iterations in a loop and then multiply the edge
weights by that number. Moreover, the workload of each thread is also
implicitly taken into account if the loop scheduling is static. This is done
by dividing the number of iterations by the number of threads. For instance,
if a loop has 100 iterations, and it is statically scheduled among four
threads, we roughly assume each thread executes 25 iterations; therefore,
edges within the body of the loop will have a weight of 25.

  * –

If statements: The nodes and edges under one branch of an if statement versus
the other will not have the same number of executions. Therefore, there is a
need for justification of edges of if statements. To illustrate the execution
of branches, we posit a probability. The assumption is that each branch of an
if statement has a probability of $\frac{1}{2}$. Therefore, the weights of
edges within if statements are divided by 2. Experimental results show that
this representation of if statements provides better results, further
optimization of edge weights is left as a future work.

###  III-B Benefit to GNN

$ParaGraph$ can be used by existing GNN models, which are a type of neural
network that can operate over graphs, for downstream tasks. Especially, those
tasks that are related to performance optimizations. $ParaGraph$ can be
exploited both as a heterogeneous graph or as a homogeneous graph. For
simplicity, in this study, we treat $ParaGraph$ as a homogeneous graph.
Typically a graph is shown as Equation 1, where $V$ is a set of nodes and $E$
contains an adjacency matrix. Elements of $E$ are either 1 or 0 to show
whether an edge exists between two nodes.

$G=(V,E)$ (1)

We extend the AST for our new representation by including new edges like
NextToken, NextSib, and other types that were explained earlier. We formally
define $ParaGraph$ as follows:

$ParaGraph=(V,E,T,W)$ (2)

Where $V$ and $E$ are previously defined in Equation 1 and
$T\in\mathbb{Z}^{+}$ represents the type of the edges (such as NextToken,
NextSib etc). $W\in\mathbb{Z}^{+}$ presents the weight which is zero for any
edge type other than Child.

These additional edges and attributes are all added to the graphs statically.
Therefore one obvious question could be why we do not let the model reason
about these additional attributes by itself rather than presenting them in the
graph? On one hand, adding these attributes to the graphs is not costly; on
the other hand, we do not want to waste the resources and the learning
capabilities of the model to realize these obvious facts. Therefore by
augmenting the ASTs and including the attributes mentioned earlier, we can
train GNN models more effectively. Later on, our ablation study further
improves that the model benefits from the attributes and information presented
in $ParaGraph$.

We adapt Relational Graph Attention Networks (RGAT) [27] and use $ParaGraph$
representation as input to train a model for predicting runtime of
applications on different accelerators. In RGAT, attention logits are computed
per each edge type. In the result section, we will see that with $ParaGraph$,
the model has a very small error in its predictions and outperforms the
current state-of-the-art approach.

Figure 3 shows the overall workflow of our GNN-based pipeline for runtime
prediction. The first step is preparing different variants of an application
for each accelerator used in this study. This step is further explained in
detail in the next section. Then, using Clang, ASTs are produced and a series
of augmentations, as discussed, are applied to them to construct $ParaGraph$.
In order to train a model to predict the runtime of an application, we need to
create a dataset with a list of applications and their variant accompanied by
their runtime. Therefore, we execute each one of the variants on the specified
accelerator and measure the runtime. Lastly, the $ParaGraph$ representation of
variants and their runtime are used to train the GNN model. Along with
$ParaGraph$s, our feature set also includes the number of teams and threads
used for executing an application. These features are then fed into the GNN-
based model for predicting the runtime.

##  IV Experiments and Setups

We used two clusters and compilers to test and evaluate our tool. Our
experiments were carried out on the ORNL’s Summit supercomputing cluster [28]
with LLVM/Clang (ver 13.0) and a nvptx backend for GPUs, as well as the LLNL’s
Corona cluster [29] with LLVM/Clang (ver 15.0) and a rocm backend for GPUs.
For the purposes of this study, we only use one GPU per cluster node.

One of the primary obstacles we faced, as with any other data-driven approach,
is the lack of a publicly accessible dataset. Creating a dataset that can be
used to train our model was the first step in building a data-driven based
cost model. Although there are plenty of OpenMP benchmark suites publicly
available, little effort has been put into OpenMP GPU offloading benchmarking.
There are very few publicly available benchmarks for OpenMP GPU offloading.
Even though a few benchmarks, such as the Rodinia benchmark suite [30],
include kernels that implement OpenMP GPU offloading, we cannot use them due
to a lack of different kernel variations. Therefore, as a first step, we
create a collection of benchmarks that leverage OpenMP GPU offloading. The
goal is to include a broad class of benchmarks that would cover a spectrum of
statistical simulation, linear algebra, data mining, etc.

###  IV-A Data Collection

There are three parts to our data collection: Code Variant Generation, Graph
Generation, and Runtime Collection.

Application |  Num Kernels |  Domain
---|---|---
Correlation Coefficient [31] |  1 |  Statistics
Covariance |  2 |  Probability Theory
Gauss Seidel Method |  1 |  Linear Algebra
K-nearest neighbors [30] |  1 |  Data Mining
Laplace’s Equation [32] |  2 |  Numerical Analysis
Matrix-Matrix Multiplication |  1 |  Linear Algebra
Matrix-Vector Multiplication |  1 |  Linear Algebra
Matrix Transpose |  1 |  Linear Algebra
Particle Filter [30] |  7 |  Medical Imaging
TABLE I: Benchmark Applications

####  IV-A1 Code Variant Generation

The first step is to create different kernels for which data needs to be
collected. We needed application kernels from multiple domains to have a
better prediction across diverse disciplines of Computer Science. Table I
specifies nine benchmark applications that we have identified for data
collection. However, nine applications are not enough. There are only
seventeen kernels if we count the number of kernels in each of these
applications. In the future, as this model develops, more applications will be
added.

As a proof of concept for this research, we only consider the following six
transformations:

  * •

cpu: A cpu parallel kernel using omp parallel for.

  * •

cpu_collapse: In case of a nested collapsible loop, we collapse it with omp
parallel for collapse(2) directive.

  * •

gpu: A gpu kernel using a combined omp target teams distribute parallel for
directive. All data is already considered to be present on the GPU.

  * •

gpu_collapse: A gpu kernel with nested collapsible loop using omp target teams
distribute parallel for collapse(2) directive. All data is already considered
to be present on the GPU.

  * •

gpu_mem: Same as combined gpu offloading, but with data transfer.

  * •

gpu_collapse_mem: Same as combined gpu_collapse, but with data transfer.

We used the OpenMP Advisor tool [20] to generate these kernel variants. Even
after creating multiple variations for each kernel, we only had 66 distinct
kernels. To create more kernels, we vary the levels of parallelism and data
used for all these kernels. We modified these parameters for each application
to generate a minimum of 2000-3000 distinct kernels for our dataset. When the
kernels from all applications are combined, we get almost 26,000 unique pieces
of information for our dataset. This step is independent of the architecture
on which we will be running the kernels.

####  IV-A2 ParaGraph Generation

In this step of our pipeline, we create the $ParaGraphs$ for the kernels
generated in the previous section. For each kernel that we want to train, our
dataset should include two parts – the GNN graph and the runtime. For
training, the first step is to generate $ParaGraph$ for each one of the
kernels. Since $ParaGraph$ is based on AST, each kernel is compiled, and ASTs
are generated. In order to produce $ParaGraph$, each AST is traversed, and
additional edges and weights are added as explained in Section III. While edge
type conveys the information from the compiler’s point of view, edge weights
expose information about the flow of execution and, more importantly, the
hotspots of the kernels, such as loops. Typically, edges within a ForLoop have
more weights to represent the number of iterations of a loop.

####  IV-A3 Runtime Collection

As our next step, we collect the actual runtime of each kernel. This process
is entirely dependent on the hardware that the kernel is intended to run
on.The complexity of this process was determined by our opportunity to get
access to the compute nodes of the HPC clusters. Failures in this step were
difficult to handle interactively, and we had to wait for individual jobs to
fail before even realizing the error. As a result, we had to be very cautious
when submitting jobs to run on these clusters.

Platform |  #Data Points |  Runtime Range (ms) |  Std. Dev.
---|---|---|---
Summit |  |  |
IBM POWER9 (CPU) |  13,023 |  [0.23 - 736,798] |  48,502
NVIDIA V100 (GPU) |  26,040 |  [0.035 - 30,174] |  3,708
Corona |  |  |
AMD EPYC7401 (CPU) |  17,681 |  [0.024 - 291,627] |  16,942
AMD MI50 (GPU) |  26,668 |  [0.448 - 46,913] |  4,828
TABLE II: Data points collected on each accelerator.

Using the OpenMP Advisor, we were able to generate code that already collected
kernel runtime information. It accomplished this by including the gettimeofday
function call, which gets the clock time of a system expressed in elapsed
seconds and microseconds, both before and after a kernel call, and then
computes the difference between the times to obtain the runtime in
microseconds. OpenMP Advisor also aided us in generating all six variants of
each kernel, and we were able to generate multiple variants for each
application by varying the levels of parallelism and data used. We built each
of these kernels on a local cluster to ensure that our code does not break as
a result of using an automated tool to generate our kernels. After
successfully building all of our kernels, we ran each kernel on a local
cluster to validate our data collection. Once we were confident that all of
our kernels were building and running on a local cluster, we created jobs to
build and run them on the Summit and Corona clusters. Nonetheless, we had
limited access to these clusters, and our job would not run for long due to
node failure or time constraints. As a result, we had to submit multiple jobs
and keep track of their success or failure in collecting data.

As soon as we had each kernel’s runtime, we could easily accompany them with
the corresponding $ParaGraphs$ to collect data points for our dataset. Table
II shows the statistics of data points collected from each one of the
accelerators.

###  IV-B ParaGraph Model

Once we have enough data on both the Summit and Corona clusters, we begin
training our model. We implemented a GNN-based neural network using RGAT [27]
as the convolution layers. Our model is implemented using Pytorch-Geometric
library with Mean Squared Error as the loss function, and Adam [33] as the
optimizer. To embed the graph, the model uses three graph convolution layers
based on RGAT, followed by two fully connected layers with ReLu activation
function. As mentioned, the number of teams and threads are considered as two
additional features. Another fully connected layer is used to embed these two
features. Finally, the embedding of the graph and the two features are
concatenated together and passed through the last fully connected layer for
runtime prediction.

The edge weights and two additional features are normalized using
MinMaxScaler. The dataset is split into train-validation sets using a 9:1
ratio.

##  V Results

In the following subsection, we explain the metrics we have used to evaluate
the $ParaGraph$ model.

###  V-A Evaluation Metric

To evaluate the performance of $ParaGraph$, we use $RMSE$ which is Root Mean
Square Error (Equation 3).

$RMSE=\sqrt{\frac{\sum_{i=1}^{N}(x_{i}-\hat{x}_{i})^{2}}{N}}$ (3)

Where $x_{i}$ stands for the runtime of a data point in microseconds,
$\hat{x}_{i}$ is the predicted runtime by $ParaGraph$, and $N$ is the total
number of samples.

Since the range of runtime differs across platforms, the normalized version of
$RMSE$ is also considered. Normalized $RMSE$ is calculated by dividing the
$RMSE$ by the distance between the minimum and maximum runtime. We also use
relative error (i.e., absolute error divided by the range of runtime) to
report the error rate.

Platform | RMSE (ms) | Norm-RMSE
---|---|---
Summit |  |
IBM POWER9 (CPU) | $4325$ | $6\times 10^{-3}$
NVIDIA V100 (GPU) | $280$ | $9\times 10^{-3}$
Corona |  |
AMD EPYC7401 (CPU) | $968$ | $4\times 10^{-3}$
AMD MI50 (GPU) | $510$ | $1\times 10^{-2}$
TABLE III: Experimental Results Figure 4: Prediction error per 10-second bins.

###  V-B Experimental Results

Table III shows the experimental results for each accelerator. We have NVIDIA
V100, AMD MI50, and IBM POWER9 with 22 cores and AMD EPYC 7401 with 24 cores.
As shown, the $RMSE$ values range from 280 (ms) to 4325 (ms). The reason why
we have different $RMSE$ values, such as 4325 (ms) for POWER9 accelerator,
lies in the fact that the runtime dispersion differs across the accelerators.
Standard Deviation in Table II gives us some insights on how dispersed the
collected data are. Moving on to Normalized-RMSE, which is independent of the
range of runtime, we see $ParaGraph$ has relatively the same error across
accelerators; thus, it can be applied to different accelerators.

To further analyze our results, we have calculated a relative error per bins
of 10 seconds. Figure 4 shows 11 bins for all four accelerators. Each bin has
a 10-second range except the last bin. The figure shows that the relative
error is small across different bins and accelerators (less than %10), meaning
that our $ParaGraph$ model has stable behavior for different ranges of
runtimes.

Moreover, we analyze how stable the model is during the training process in
Figure 5. The figure shows validation $RMSE$ for four accelerators. In the
first few epochs, the $ParaGraph$ model is not very stable, resulting in
fluctuations in the $RMSE$; however, as the model is trained further, it is
able to pick and learn the features from the representation and reduce $RMSE$
value per each epoch and ultimately converges.

Figure 5: Normalized RMSE per each epoch. Figure 6: Error rate per each
application.

Lastly, we calculate the average relative error per application to see if the
$ParaGraph$ model is able to have prediction with low error for all types of
applications. Figure 6 shows the error rate of each application. As can be
seen, the model is indeed capable of having accurate predictions for a wide
variety of applications resulting in a low error rate. Therefore, the model is
not biased toward one application. On the AMD MI50 GPU, the Laplace data was
corrupted during collection. Consequently, neither this study nor the training
process includes that data.

###  V-C Ablation Study

$ParaGraph$ representation, as explained in Section III, is built by applying
a series of major augmentations on top of AST. In this section, we quantify
the impact of these augmentations. First, we consider the AST itself without
additional edges and weights, we call it Raw AST. Then, we add additional
edges and edge types to the AST and name it Augmented AST. Lastly, we have
$ParaGraph$ that contains both additional edges and also edge weights. Table
IV shows the results of the ablation study. We see that Raw AST results in the
highest error for all four accelerators. Adding new edges and introducing new
types of edges (Augmented AST) improve the prediction to some extent. For
example, the $RMSE$ of V100 drops from 2114 (ms) to 786 (ms) with the addition
of these new edges.

One of the key characteristics of our proposed program representation is the
edge weights. Edge weights convey essential information about how often
different regions of the AST execute. Therefore, we see quite a good
improvement in $RMSE$ when edge weights are added. For instance, the $RMSE$
for V100 is further improved to 280 (ms).

Platform |  Raw AST |  Aug AST |  ParaGraph
---|---|---|---
Summit |  |  |
IBM POWER9 (CPU) |  $27593$ |  $26860$ |  4325
NVIDIA V100 (GPU) |  2114 |  786 |  280
Corona |  |  |
AMD EPYC7401 (CPU) |  $11911$ |  $9633$ |  968
AMD MI50 (GPU) |  $2888$ |  $1177$ |  510
TABLE IV: RMSE of training with and without edges’ weight. Figure 7: RMSE of
the validation set during training the GNN models on MI50 data points.

We analyze the addition of edges and their weights further to see how the
training process of the $ParaGraph$ model is affected by these augmentations.

Figure 7 shows how the model behaves during training. This figure depicts the
$RMSE$ value per each epoch for Raw AST, Augmented AST, and $ParaGraph$ on the
MI50 accelerator.

Using only the Raw AST without any augmentations, which means having only one
edge type, the model is able to learn some characteristics of the applications
and reduce the $RMSE$ per epoch; however, this reduction in $RMSE$ is not
significant. Augmented AST contains 8 different types of edges. We see the
addition of these edges destabilized the training process of the model. In the
first few epochs, the model is challenged to learn different relations between
the nodes. However, eventually, after several epochs, the prediction of the
model is stabilized and it achieves $RMSE$ of 1177 (ms). Once edges of the
Augmented AST are augmented with weights, thus $ParaGraph$ is constructed, we
see further improvements in the model’s predictions. Although the validation
$RMSE$ has some fluctuation in the initial epochs, it ultimately converges
with a considerably smaller error.

###  V-D Comparison with State-of-the-art Tool

To the best of our knowledge, $COMPOFF$ [2] is the only state-of-the-art
OpenMP GPU offloading cost model. We further compare the results from
$ParaGraph$ with those of $COMPOFF$. As mentioned in Section II-D OpenMP
Advisor uses $COMPOFF$ for predicting runtime for OpenMP kernels.

While OpenMP Advisor eventually needs a cost model that can forecast for all
possible underlying architectures, $COMPOFF$ is currently only suitable for
GPU execution. In contrast, $ParaGraph$ can model the application runtime
regardless of the underlying architecture. In our experiments, $ParaGraph$ was
used to predict runtime on both CPUs and GPUs. Therefore, we compare
$ParaGraph$ and $COMPOFF$ using data points collected from the NVIDIA V100 GPU
only.

For both $ParaGraph$ and $COMPOFF$, Figure 9 demonstrates a strong correlation
between the actual and predicted data. The results for $COMPOFF$ are
represented by blue dots, while those for $ParaGraph$ are represented by
orange dots. Despite the fact that both perform admirably, $ParaGraph$
demonstrates a much stronger correlation between the predicted and actual
runtime. In Figure 8, $COMPOFF$ (blue) demonstrates a slightly higher error
rate for smaller runtime kernels, but as the runtime increases, this error
rate decreases (or disappears entirely). However, for all kernels, the error
rate is significantly lower for $ParaGraph$ (orange). Similar results were
also observed on the AMD GPU. We were unable to compare the outcomes of our
prediction on CPUs because $COMPOFF$ does not support CPUs. This is one
significant advantage $ParaGraph$ has over $COMPOFF$.

Figure 8: Comparison of ParaGraph and COMPOFF on predicting runtimes on NVIDIA
V100 for each data point.

###  V-E Discussion

$ParaGraph$ is able to model the execution of applications statically. Our
experimental results show the effectiveness of $ParaGraph$ in predicting
applications’ runtime across different accelerators. On the other hand, there
are applications whose runtime behavior is not stable. That is, they show
different behaviors across various executions. These applications can not be
modeled statically. As a result, predicting their runtime using only static
features is a challenge. One way of tackling this challenge could be
developing a hybrid approach using static and dynamic features (for example,
Performance Counters) of applications. However, collecting dynamic features of
applications requires the execution of applications which will potentially
increase the cost of prediction.

Figure 9: ParaGraph and COMPOFF prediction on NVIDIA V100 as compared to the
actual runtime.

##  VI Conclusion and Future Work

In this paper, we proposed $ParaGraph$. A novel way of representing HPC
kernels. $ParaGraph$ is based on ASTs; it incorporates additional features to
the AST to better represent applications. In particular, $ParaGraph$ includes
some of the insights of the compiler by adding new edges to the AST.
$ParaGraph$ also adds information about the way loops and if statements will
be executed by adding weights to the corresponding edges. We evaluated
$ParaGraph$ on a set of applications for four different accelerators. It
achieved a low error rate highlighting the effectiveness of $ParaGraph$. One
of the main benefits of $ParaGraph$ is its independence of underlying hardware
architecture. Therefore, it can be applied to a wider range of architectures,
unlike previous approaches.

In this work, we used $ParaGraph$ for predicting the execution time of a given
kernel for OpenMP GPU offloading problem. We plan to explore and analyze how
$ParaGraph$ can help other OpenMP optimization strategies such as predicting
SIMD stride, scheduling strategies, loop chunk size, etc. in the future.
Another interesting research area to explore is to use $ParaGraph$ and capture
parallelism for other parallel programming models such as OpenACC, Kokkos,
HIP, etc. for building cost models for various optimization problems.

## References

  * [1] L. Dagum and R. Menon, “Openmp: an industry standard api for shared-memory programming,” _IEEE computational science and engineering_ , vol. 5, no. 1, pp. 46–55, 1998.
  * [2] A. Mishra, S. Chheda, C. Soto, A. M. Malik, M. Lin, and B. Chapman, “COMPOFF: A compiler cost model using machine learning to predict the cost of openmp offloading,” in _2022 IEEE International Parallel and Distributed Processing Symposium Workshops (IPDPSW), May 30-June 3, 2022_. IEEE, 2022.
  * [3] A. Ciurumelea, S. Proksch, and H. C. Gall, “Suggesting comment completions for python using neural language models,” in _2020 IEEE 27th International Conference on Software Analysis, Evolution and Reengineering (SANER)_. IEEE, 2020, pp. 456–467.
  * [4] C. Cummins, Z. V. Fisches, T. Ben-Nun, T. Hoefler, M. F. O’Boyle, and H. Leather, “Programl: A graph-based program representation for data flow analysis and compiler optimizations,” in _International Conference on Machine Learning_. PMLR, 2021, pp. 2244–2253.
  * [5] V. Raychev, M. Vechev, and E. Yahav, “Code completion with statistical language models,” in _Proceedings of the 35th ACM SIGPLAN conference on programming language design and implementation_ , 2014, pp. 419–428.
  * [6] M. Allamanis, H. Peng, and C. Sutton, “A convolutional attention network for extreme summarization of source code,” in _International conference on machine learning_. PMLR, 2016, pp. 2091–2100.
  * [7] M. Allamanis, M. Brockschmidt, and M. Khademi, “Learning to represent programs with graphs,” _arXiv preprint arXiv:1711.00740_ , 2017.
  * [8] M. A. Heroux, R. Thakur, L. McInnes, J. S. Vetter, X. S. Li, J. Aherns, T. Munson, and K. Mohror, “Ecp software technology capability assessment report,” Oak Ridge National Lab.(ORNL), Oak Ridge, TN (United States), Tech. Rep., 2020.
  * [9] Y. Kim, J. Lee, J.-S. Kim, H. Jei, and H. Roh, “Comprehensive techniques of multi-gpu memory optimization for deep learning acceleration,” _Cluster Computing_ , vol. 23, no. 3, pp. 2193–2204, 2020.
  * [10] K. Rzadca, P. Findeisen, J. Swiderski, P. Zych, P. Broniek, J. Kusmierek, P. Nowak, B. Strack, P. Witusowski, S. Hand _et al._ , “Autopilot: workload autoscaling at google,” in _Proceedings of the Fifteenth European Conference on Computer Systems_ , 2020, pp. 1–16.
  * [11] M. Mirka, G. Sassatelli, and A. Gamatié, “Online learning for dynamic control of openmp workloads,” in _2020 9th International Conference on Modern Circuits and Systems Technologies (MOCAST)_. IEEE, 2020, pp. 1–6.
  * [12] N. Denoyelle, B. Goglin, E. Jeannot, and T. Ropars, “Data and thread placement in numa architectures: A statistical learning approach,” in _Proceedings of the 48th International Conference on Parallel Processing_ , 2019, pp. 1–10.
  * [13] J. Alcaraz, A. TehraniJamsaz, A. Dutta, A. Sikora, A. Jannesari, J. Sorribes, and E. Cesar, “Predicting number of threads using balanced datasets for openmp regions,” _Computing_ , pp. 1–19, 2022.
  * [14] A. Dutta, J. Alcaraz, A. TehraniJamsaz, A. Sikora, E. Cesar, and A. Jannesari, “Pattern-based autotuning of openmp loops using graph neural networks,” in _2022 IEEE/ACM International Workshop on Artificial Intelligence and Machine Learning for Scientific Applications (AI4S)_. IEEE, 2022, pp. 26–31.
  * [15] A. M. Malik, “Automatic static feature generation for compiler optimization problems,” in _Australasian Joint Conference on Artificial Intelligence_. Springer, 2011, pp. 769–778.
  * [16] S. Dublish, V. Nagarajan, and N. Topham, “Poise: Balancing thread-level parallelism and memory system performance in gpus using machine learning,” in _2019 IEEE International Symposium on High Performance Computer Architecture (HPCA)_. IEEE, 2019, pp. 492–505.
  * [17] A. Iranfar, W. S. De Souza, M. Zapater, K. Olcoz, S. X. de Souza, and D. Atienza, “A machine learning-based framework for throughput estimation of time-varying applications in multi-core servers,” in _2019 IFIP/IEEE 27th International Conference on Very Large Scale Integration (VLSI-SoC)_. IEEE, 2019, pp. 211–216.
  * [18] D. Grewe and M. F. O’Boyle, “A static task partitioning approach for heterogeneous systems using opencl,” in _International conference on compiler construction_. Springer, 2011, pp. 286–305.
  * [19] H. Sayadi, “Energy-efficiency prediction of multithreaded workloads on heterogeneous composite cores architectures using machine learning techniques,” _arXiv preprint arXiv:1808.01728_ , 2018.
  * [20] A. Mishra, A. M. Malik, M. Lin, and B. Chapman, “Openmp advisor,” _arXiv preprint arXiv:2301.03636_ , 2023.
  * [21] R. B. Roy, T. Patel, V. Gadepally, and D. Tiwari, “Bliss: auto-tuning complex applications using a pool of diverse lightweight learning models,” in _Proceedings of the 42nd ACM SIGPLAN International Conference on Programming Language Design and Implementation_ , 2021, pp. 1280–1295.
  * [22] S. Falkner, A. Klein, and F. Hutter, “Bohb: Robust and efficient hyperparameter optimization at scale,” in _International Conference on Machine Learning_. PMLR, 2018, pp. 1437–1446.
  * [23] H. Menon, A. Bhatele, and T. Gamblin, “Auto-tuning parameter choices in hpc applications using bayesian optimization,” in _2020 IEEE International Parallel and Distributed Processing Symposium (IPDPS)_. IEEE, 2020, pp. 831–840.
  * [24] Y. Liu, W. M. Sid-Lakhdar, O. Marques, X. Zhu, C. Meng, J. W. Demmel, and X. S. Li, “Gptune: multitask learning for autotuning exascale applications,” in _Proceedings of the 26th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming_ , 2021, pp. 234–246.
  * [25] X. Wu, M. Kruse, P. Balaprakash, H. Finkel, P. Hovland, V. Taylor, and M. Hall, “Autotuning polybench benchmarks with llvm clang/polly loop optimization pragmas using bayesian optimization,” _Concurrency and Computation: Practice and Experience_ , vol. 34, no. 20, p. e6683, 2022.
  * [26] M. Allamanis, “Graph neural networks in program analysis,” in _Graph Neural Networks: Foundations, Frontiers, and Applications_. Springer, 2022, pp. 483–497.
  * [27] D. Busbridge, D. Sherburn, P. Cavallo, and N. Y. Hammerla, “Relational graph attention networks,” _arXiv preprint arXiv:1904.05811_ , 2019.
  * [28] ORNL, “Oak Ridge Leadership Computing Facility - Summit supercomputing cluster,” 2017. [Online]. Available: https://www.olcf.ornl.gov/summit/
  * [29] L. L. N. Laboratory, “LLNL - Corona,” 2019. [Online]. Available: https://hpc.llnl.gov/hardware/compute-platforms/corona
  * [30] S. Che, M. Boyer, J. Meng, D. Tarjan, J. W. Sheaffer, S.-H. Lee, and K. Skadron, “Rodinia: A benchmark suite for heterogeneous computing,” in _2009 IEEE international symposium on workload characterization (IISWC)_. Ieee, 2009, pp. 44–54.
  * [31] P. Schober, C. Boer, and L. A. Schwarte, “Correlation coefficients: appropriate use and interpretation,” _Anesthesia & Analgesia_, vol. 126, no. 5, pp. 1763–1768, 2018.
  * [32] A. K. Mitra, “Finite difference method for the solution of laplace equation,” _Department of aerospace engineering Iowa state University_ , 2010.
  * [33] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.


# Thermodynamic uncertainty relations in the presence of non-linear friction
and memory

A. Plati1 A. Puglisi2,3,4 A. Sarracino5

(1Université Paris-Saclay, CNRS, Laboratoire de Physique des Solides, 91405
Orsay, France
2Department of Physics, University of Rome Sapienza, P.le Aldo Moro 2, 00185,
Rome, Italy
3Institute for Complex Systems - CNR, P.le Aldo Moro 2, 00185, Rome, Italy
4INFN, University of Rome Tor Vergata, Via della Ricerca Scientifica 1, 00133,
Rome, Italy
5Department of Engineering, University of Campania “Luigi Vanvitelli”, 81031
Aversa (CE), Italy
)

###### Abstract

A new Thermodynamic Uncertainty Relation (TUR) is derived for systems
described by linearly coupled Langevin equations in the presence of non-linear
frictional forces. In our scheme, the main variable represents the velocity of
a particle, while the other coupled variables describe memory effects which
may arise from strongly correlated degrees of freedom with several time-scales
and, in general, are associated with thermal baths at different temperatures.
The new TUR gives a lower bound for the mean-squared displacement of the
position of the particle, including its asymptotic diffusion coefficient. This
bound, in several examples worked out here, appears to be a good analytical
estimate of the real diffusion coefficient. The new TUR can be also applied in
the absence of any external force (with or without thermal equilibrium between
the baths), a case which usually goes beyond the scope of original TURs. We
show applications to non-linear frictional models with memory, such as the
Coulomb and the Prantdtl-Tomlinson models, usually representative of friction
at the nano-scale and within atomic-force microscopy experiments.

##  1 Introduction

In the last decade a family of rigorous inequalities has been derived,
generally known as Thermodynamic Uncertainty Relations (TURs) [12]. They
provide bounds to fluctuations of several quantities, and those bounds are
usually related to entropy production. The most simple and used result can be
summarised as the entropy production rate limiting from above the precision
rate of any non-equilibrium current of a driven system (a quantity similar but
not exactly identical to a signal-to-noise ratio) [1, 8, 25, 27, 10]. The TURs
constitute a rare general result for out-of-equilibrium statistical mechanics,
akin to Fluctuation-Dissipation Relations or Fluctuation Relations [18]. They
also provide an estimate of quantities which are not easily measurable, and
this occurs in a twofold way: in certain cases entropy production is
accessible/predictable, while the currents are not; in other cases the
opposite happens [9]. Of course, none of those problems is immediately solved
by an inequality, but several examples have been given where the inequality is
not far from being saturated [13]. Moreover, a bound can be exploited through
optimisation strategies. For instance, entropy production has been
successfully approximated by machine-learning-based protocols, looking for the
most precise non-equilibrium current which approaches the bound [14].

In this paper we build upon a recently proposed inequality that relates
entropy production and the mean-squared displacement (MSD) of a diffusing
particle, valid also when in the presence of memory effects (e.g. degrees of
freedom which are coupled to the particle velocity with different time-scales)
[20]. The main novelty in this work is represented by the possibility of
considering also non-linear friction forces, i.e. dissipative forces that
depend also non-linearly on the particle velocity [7, 11, 21, 2, 24, 17, 4].
This category includes typical models such as the classical Coulomb friction,
as well as power laws observed in foams [3] and other non-linear models (such
as the rate-and-state equation) that occur when considering microscopic and
nanoscopic friction experiments [28]. We will show how, in this framework, it
is possible to derive analytically a tight lower bound for the diffusion
coefficient, a quantity that requires the analytical form of the propagator to
be computed exactly [15].

The paper is organized as follows. In Section 2 we discuss a general strategy
to derive TURs which is different from previously proposed approaches, and
which is particularly useful in our case. In Section 3, we apply such a
derivation to the problem of generalised Langevin equations with exponential
memory kernel and non-linear frictional velocity-dependent forces. This class
of systems can be mapped - by Markovian embedding - into linearly coupled
Markovian Langevin equations, still maintaining the non-linear forces in the
equation for the particle motion. In Section 4, we draw conclusions and
propose future perspectives.

##  2 General strategies for the derivation of TURs

The TURs have been initially conjectured and then derived, by means of several
techniques, for a larger and larger category of systems, encompassing
Markovian discrete and continuous dynamics as well as transient and steady
states [1, 8, 6]. In recent years, a general strategy to derive TURs-like
inequalities, including bounds for linear response, has been traced in the
framework of information theory [5]. It goes through the introduction of a
virtual perturbation in the original model equations, and then the application
of the Cramér-Rao inequality, which identifies the reciprocal of the Fisher
Information as a lower bound for the variance of an unbiased estimator [26].
The usual application of this idea is introducing the virtual perturbation in
the equations with the aim of modifying the irreversible probability current
of the system, that is the part of the probability current which is
responsible for the lack of detailed balance (and is therefore connected to
entropy production). However, such a strategy cannot be directly applied to a
system with underdamped dynamics, raising the necessity of alternative
treatments [10]. Recently, we have proposed a different application of the
Cramér-Rao-based derivation, in the context of Stochastic Differential
Equations (SDE): the virtual perturbation is applied directly to the variables
appearing in the SDE, and the symmetries of the system are exploited to have
additional contributions that can be optimised to obtain tighter bounds [20].
Here we build upon such an approach to include a non-linear friction-force
term in the main (velocity) variable.

###  2.1 General bound in the presence of non-linear friction and memory

We consider a Generalized Langevin Equation with exponential memory kernel
[16] for a particle with position $\theta$ and velocity $\omega$, subjected to
a nonlinear force $\mathcal{F}(\omega)$ and dragged by a constant force
$F_{\text{ext}}$ [22]:

$\displaystyle\dot{\theta}$ $\displaystyle=$ $\displaystyle\omega,$ (1)
$\displaystyle\dot{\omega}(t)$ $\displaystyle=$
$\displaystyle-\int_{-\infty}^{t}\gamma(t-t^{\prime})\omega(t^{\prime})dt^{\prime}+\eta_{s}(t)+F_{\text{ext}}+\mathcal{F}(\omega),$
(2) $\displaystyle\gamma(t)$ $\displaystyle=$
$\displaystyle\frac{2}{\tau}\delta(t)+\sum_{k=1}^{n}a_{k}e^{-\frac{t-t^{\prime}}{\tau_{k}}},$
(3) $\displaystyle\langle\eta_{s}(t)\eta_{s}(t^{\prime})\rangle$
$\displaystyle=$
$\displaystyle\frac{2q}{\tau}\delta(|t-t^{\prime}|)+\sum_{k=1}^{n}q_{k}a_{k}e^{-\frac{|t-t^{\prime}|}{\tau_{k}}}.$
(4)

Upon defining the $n$ auxiliary variables:

$\Omega_{k}=-b_{k}\int_{-\infty}^{t}dt^{\prime}e^{-\frac{t-t^{\prime}}{\tau_{k}}}\left[a_{k}\omega(t^{\prime})-\sqrt{\frac{2q_{k}a_{k}}{\tau_{k}}}\xi_{k}(t^{\prime})\right]\quad
k\in[1,n],$ (5)

the system can be written as a Markovian SDE for the $n+1$ dimensional vector
$\boldsymbol{X}=\\{\omega,\Omega_{1},\ldots,\Omega_{n}\\}$:

$\dot{\omega}=-\frac{\omega}{\tau}+\sum_{k}\frac{\Omega_{k}}{b_{k}}+\sqrt{\frac{2q}{\tau}}\xi(t)+F_{\text{ext}}+\mathcal{F}(\omega),$
(6a)
$\dot{\Omega}_{k}=-\frac{\Omega_{k}}{\tau_{k}}-a_{k}b_{k}\omega+\sqrt{\frac{2q_{k}a_{k}b_{k}^{2}}{\tau_{k}}}\xi_{k}(t).$
(6b)

We define accordingly the noise amplitudes $B_{i}$, such as
$B_{0}=\sqrt{2q/\tau}$ and $B_{i}=\sqrt{2q_{i}a_{i}b_{i}^{2}/\tau_{i}}$ for
$i>0$. From Eqs. (6) we obtain the following relation between the average
values of $\omega$ and $\Omega_{k}$s in the stationary state:

$\quad\langle\Omega_{k}\rangle=-\tau_{k}a_{k}b_{k}\langle\omega\rangle.$ (7)

The strategy to derive the new inequality is to introduce a perturbation in
Eqs. (6) by adding to the r.h.s. of each equation the corresponding component
of the vector $h\boldsymbol{V}(\omega)$, which is defined as:

$V_{0}(\omega)=\langle\omega\rangle/\tau-\langle\omega\rangle\mathcal{F}^{\prime}(\omega)-\sum_{k}u_{k}/b_{k}\quad\text{and}\quad
V_{k}=a_{k}b_{k}\langle\omega\rangle+u_{k}/\tau_{k},$ (8)

where $u_{k}$ are arbitrary parameters, to be fixed later.

Considering that the $h$-perturbation can be arbitrarily small, we are allowed
to take
$\mathcal{F}(\omega)-h\langle\omega\rangle\mathcal{F}^{\prime}(\omega)\simeq\mathcal{F}(\omega-h\langle\omega\rangle)$.
Then, the perturbed dynamics is equivalent to the dynamics described by Eqs.
(6) for the vector
$\boldsymbol{X}_{h}=\\{\omega-h\langle\omega\rangle,\Omega_{1}-hu_{1},\ldots,\Omega_{n}-hu_{n}\\}$.
In other words, the perturbed $h-$dynamics is the same as the original one,
but with shifted averaged values:

$\langle\omega\rangle_{h}=(1+h)\langle\omega\rangle,\quad\langle\Omega_{k}\rangle_{h}=\langle\Omega_{k}\rangle+hu_{k},$
(9)

and for the perturbed steady state we have
$P_{h}(\boldsymbol{X})=P(\boldsymbol{X}_{h})$, where $P_{h}$ and $P$ are the
stationary distributions for the perturbed and unperturbed dynamics,
respectively. Now we apply the Cramér-Rao inequality that, in our case, reads:

$\frac{\text{Var}_{h=0}(\theta(t))}{[\partial_{h=0}\langle\theta(t)\rangle_{h}]^{2}}=\frac{\langle\Delta\theta(t)^{2}\rangle}{(\langle\omega\rangle
t)^{2}}\geq\frac{1}{\mathcal{I}_{\text{F}}(h=0)},$ (10)

where $\Delta\theta(t)=\theta(t)-\langle\theta\rangle$,

$\mathcal{I_{\text{F}}}(h=0)=-\langle\partial^{2}_{h}\ln
P_{h}(\boldsymbol{X})\rangle_{h=0}+\left\langle\int_{0}^{t}dt^{\prime}\sum_{i}\left(\frac{V_{i}(\omega)}{B_{i}}\right)^{2}\right\rangle_{h=0}=\mathcal{I}/2+\mathcal{B}t,$
(11)

and

$\mathcal{I}=2\int
d\boldsymbol{X}\frac{[\partial_{h}P_{h}(\boldsymbol{X})]^{2}|_{h=0}}{P(\boldsymbol{X})},\quad\mathcal{B}=\left(\frac{\langle\omega\rangle}{\tau}-\langle\omega\rangle\langle\mathcal{F}^{\prime}\rangle-\sum_{j}\frac{u_{j}}{b_{j}}\right)^{2}B_{0}^{-2}+\sum_{k}\left(a_{k}b_{k}\langle\omega\rangle+\frac{u_{k}}{\tau_{k}}\right)^{2}B_{k}^{-2}.$
(12)

Note that $\langle\mathcal{F}^{\prime}\rangle$ requires the knowledge of the
stationary distribution $P(\boldsymbol{X}_{h})$ to be computed. The values of
$u_{k}$ as a function of the model parameters can be fixed by minimising
$\mathcal{B}(\boldsymbol{u})$. In this way, one obtains a bound which is
optimal in this set of possible perturbations. This is done by solving the
system of equations defined by $\partial\mathcal{B}(\boldsymbol{u})/\partial
u_{k}=0$ $\forall k$. The system is formally solved by
$\boldsymbol{u}=\hat{T}^{-1}\boldsymbol{z}$ where:

$T_{ii}=\frac{1}{\tau_{i}^{2}B_{i}^{2}}+\frac{1}{b_{i}^{2}B_{0}^{2}},\quad
T_{i\neq j}=\frac{1}{b_{i}b_{j}B_{0}^{2}},\quad
z_{i}=\langle\omega\rangle\left[\frac{\tau^{-1}+\langle\mathcal{F}^{\prime}\rangle}{b_{i}B_{0}^{2}}-\frac{a_{i}b_{i}}{\tau_{i}B_{i}^{2}}\right].$
(13)

First, we note that with equal thermostats $q_{i}=q$ $\forall i$ and
$\langle\mathcal{F}^{\prime}\rangle=0$, one has $z_{i}=0$, which means that
the optimal solution is $u_{k}=0$ $\forall$ $k$. This is consistent with the
analytical result for the linear system where the bound for the diffusion
coefficient is saturated at equilibrium.

In order to find the analytical expression of the set of values $\\{u_{k}\\}$,
an exact inversion of the matrix $\hat{T}$ is required, leading to very
cumbersome formulae. However, from Eq. (13), one has that the solutions can be
written as $u_{k}=\langle\omega\rangle\tilde{u}_{k}$. It is important to note
that this is not an assumption of linearity in $\langle\omega\rangle$, but a
simple rescaling that will allow us to obtain some simplifications in the
final inequality.

Now we can define $\mathcal{I}=\langle\omega\rangle^{2}\tilde{\mathcal{I}}$
and $\mathcal{B}=\langle\omega\rangle^{2}\tilde{\mathcal{B}}$ and simplify
$\langle\omega\rangle^{2}$ in the bound for the msd Eq. (10) obtaining:

$\langle\Delta\theta(t)^{2}\rangle\geq\frac{t^{2}}{\tilde{\mathcal{B}}t+\tilde{\mathcal{I}}/2}.$
(14)

We note that, from Eq. (12), it is easy to verify that $\mathcal{I}$ is
proportional to $\langle\omega\rangle^{2}$. Note that from the equation above
and Eqs. (12), one has that the case without memory ($a_{k}=u_{k}=0$ $\forall$
$k$) could lead to the pathological condition
$1/\tau-\langle\mathcal{F}^{\prime}\rangle=0$ for which Eq. (14) predicts a
non-physical asymptotic ballistic regime. However, for frictional forces,
which constitute the main focus of this paper, $\mathcal{F}^{\prime}$ is non-
positive, thus preventing any irregularity of the bound. Anyway, even when
$\mathcal{F}^{\prime}\geq 0$, the system without memory admits a potential
solution so that the quantity $1/\tau-\langle\mathcal{F}^{\prime}\rangle$ can
be evaluated analytically to check the validity of Eq. (14). Moreover, from
Eqs. (12) and (14), one can also check that, in the non-optimal case
($u_{k}=0$ $\forall$ $k$) and without nonlinear friction
($\mathcal{F}(\omega)=0$), the TUR derived in Ref. [20] for linear forces is
recovered. In this condition, the specific contribution to the entropy
production rate due to the presence of external forcing appears in the
denominator of the bound.

Both $\tilde{\mathcal{B}}$ and $\tilde{\mathcal{I}}$ are not known in general,
since they require the knowledge of the stationary distribution
$P(\boldsymbol{X})$. In the following, we will concentrate on cases where the
r.h.s. of Eq. (14) can be found analytically.

###  2.2 Two classes of models

Equation (14) can be used as an analytical bound for the diffusion coefficient
in the following classes of models: i) models with no memory, featuring
generic nonlinear friction; ii) models with memory, including nonlinear
friction with a bounded derivative, in the large time limit. In the first
case, we have a gradient system for a single degree of freedom. We can then
find the stationary probability distribution and in turn compute
$\tilde{\mathcal{B}}$ and $\tilde{\mathcal{I}}$. In the second case,
exploiting the maximum of the modulus of the derivative of the non-linear
force $\max(|\mathcal{F}^{\prime}|)$, we can use the following chain of
inequalities

$\left(\frac{\langle\omega\rangle}{\tau}-\langle\omega\rangle\langle\mathcal{F}^{\prime}\rangle-\sum_{j}\frac{u_{j}}{b_{j}}\right)^{2}\leq\left(\frac{\langle\omega\rangle}{\tau}+\langle\omega\rangle|\langle\mathcal{F}^{\prime}\rangle|-\sum_{j}\frac{u_{j}}{b_{j}}\right)^{2}\leq\left(\frac{\langle\omega\rangle}{\tau}+\langle\omega\rangle\max(|\mathcal{F}^{\prime}|)-\sum_{j}\frac{u_{j}}{b_{j}}\right)^{2},$
(15)

and define

$\mathcal{B}_{\text{max}}=\left(\frac{\langle\omega\rangle}{\tau}+\langle\omega\rangle\max(|\mathcal{F}^{\prime}|)-\sum_{j}\frac{u_{j}}{b_{j}}\right)^{2}B_{0}^{-2}+\sum_{k}\left(a_{k}b_{k}\langle\omega\rangle+\frac{u_{k}}{\tau_{k}}\right)^{2}B_{k}^{-2}\geq\mathcal{B},$
(16)

to have an approximate but analytical bound for the MSD at long times (i.e.
for the diffusion coefficient):

$\lim_{t\to\infty}\frac{1}{t}\langle\Delta\theta(t)^{2}\rangle\geq\frac{1}{\tilde{\mathcal{B}}_{\text{max}}},$
(17)

where
$\tilde{\mathcal{B}}_{\text{max}}=\mathcal{B}_{\text{max}}/\langle\omega\rangle^{2}$.
This is an explicit bound with a closed analytical form. It is however weaker
than the non-analytic one in which $\langle\mathcal{F}^{\prime}\rangle$
explicitly appears. Nevertheless, it can still be improved by choosing the
optimal $u_{k}$ that minimize
$\tilde{\mathcal{B}}_{\text{max}}(\boldsymbol{u})$ (note that the chain of
inequalities (15) is written for generic $u_{k}$).

###  2.3 Comparison with the standard TUR

For a large class of stochastic processes, which include the model Eqs. (6)
under study, a bound for asymptotic diffusion can be derived through the
standard TUR involving the entropy production rate (EPR)
$\langle\dot{S}\rangle$ in the steady state [12, 20]:

$\lim_{t\to\infty}\frac{1}{t}\langle\Delta\theta(t)^{2}\rangle\geq\frac{2\langle\omega\rangle^{2}}{\langle\dot{S}\rangle}.$
(18)

Here, we compare this previous result with the bound Eq. (14) derived in this
paper. The EPR in the steady state for our model can be easily computed by
splitting the non-linear velocity-dependent force into a reversible and
irreversible part:
$\mathcal{F}(\omega)=\mathcal{F}_{\text{rev}}(\omega)+\mathcal{F}_{\text{irr}}(\omega)$,
where
$\mathcal{F}_{\text{rev}}(\omega)=1/2[\mathcal{F}(\omega)+\mathcal{F}(-\omega)]$
and
$\mathcal{F}_{\text{irr}}(\omega)=1/2[\mathcal{F}(\omega)-\mathcal{F}(-\omega)]$
. Then, using the results in [22], we obtain:

$\begin{split}\langle\dot{S}\rangle=&\frac{1}{q}\langle\omega
F_{\text{ext}}\rangle+\frac{1}{q}\langle\omega\mathcal{F}_{\text{rev}}(\omega)\rangle-\langle\omega\rangle^{2}\sum_{i}a_{i}\tau_{i}\left(\frac{1}{q}-\frac{1}{q_{i}}\right)+\sum_{i}\frac{1}{b_{i}}\left(\frac{1}{q}-\frac{1}{q_{i}}\right)\sigma_{0i}\\\
&-\frac{\tau}{q}\sum_{i}\frac{1}{b_{i}}\langle\mathcal{F}_{\text{irr}}(\omega)\Omega_{i}\rangle-\frac{\tau}{q}\langle\mathcal{F}_{\text{irr}}(\omega)F_{\text{ext}}\rangle-\frac{\tau}{q}\langle\mathcal{F}_{\text{irr}}(\omega)\mathcal{F}_{\text{rev}}(\omega)\rangle,\end{split}$
(19)

where
$\sigma_{0i}=\langle(\omega-\langle\omega\rangle)(\Omega_{i}-\langle\Omega_{i}\rangle)\rangle$.

In the cases of interest for this paper, the velocity-dependent force has only
the irreversible part ($\mathcal{F}_{\text{rev}}(\omega)=0$) and the
expression substantially simplifies:

$\langle\dot{S}\rangle=\frac{1}{q}F_{\text{ext}}\langle\omega\rangle-\langle\omega\rangle^{2}\sum_{i}a_{i}\tau_{i}\left(\frac{1}{q}-\frac{1}{q_{i}}\right)+\sum_{i}\frac{1}{b_{i}}\left(\frac{1}{q}-\frac{1}{q_{i}}\right)\sigma_{0i}-\frac{\tau}{q}\sum_{i}\frac{1}{b_{i}}\langle\mathcal{F}_{\text{irr}}(\omega)\Omega_{i}\rangle-\frac{\tau}{q}F_{\text{ext}}\langle\mathcal{F}_{\text{irr}}(\omega)\rangle,$
(20)

where $\mathcal{F}_{\text{irr}}(\omega)$ refers to any frictional force. The
case without memory is simpler:

$\langle\dot{S}\rangle=\frac{1}{q}F_{\text{ext}}\langle\omega\rangle-\frac{\tau}{q}F_{\text{ext}}\langle\mathcal{F}_{\text{irr}}(\omega)\rangle.$
(21)

As we expect from [23], the dynamics is reversible when there is a single
thermostat and $F_{\text{ext}}=0$. Note that in the presence of a spatial
potential, the system would have been out-of-equilibrium also in the case
$F_{\text{ext}}=0$ [23].

The bound predicted by the standard TUR coincides with Eq. (18) with
$\langle\dot{S}\rangle$ given by Eqs. (20) and (21). However, there are many
crucial differences between this bound and Eq. (14):

  1. 1.

The standard bound is valid only in the long-time regime.

  2. 2.

It requires the computation of non-trivial averages on the steady state, that
are more difficult to bound with respect to
$\langle\mathcal{F}^{\prime}\rangle$. Thus, making the case with memory
analytical is non-trivial.

  3. 3.

In general, its denominator does not scale with $\langle\omega^{2}\rangle$ and
therefore the case of free diffusion (i.e. without any applied force) becomes
more delicate to treat; in our result, on the contrary, the scaling with
$\langle\omega^{2}\rangle$ allows the free diffusion limit to have a well
defined bound.

From the above remarks, it follows that a connection with our bound is not a
direct task for the case with $F_{ext}=0$. A comparison can be done (see
Section 3.3) in the absence of memory: exploiting the potential form of the
solution of the gradient system, one can explicitly compute the averages and
check if a good scaling of $\langle\mathcal{F}_{\text{irr}}(\omega)\rangle$
with $\langle\omega\rangle$ occurs. In general, the bound derived from our
procedure is more appropriate for the estimate of the diffusion coefficient.

##  3 Applications

In this Section, we apply our results to specific models. In particular, we
consider two relevant kinds of non-linear velocity-dependent forces
$\mathcal{F}(\omega)$: the Coulomb and the Prandtl-Tomlinson (PT) friction
models. For each of them, we study a case without memory and a case with
memory with a single characteristic time. The class of considered systems is
defined by the general model:

$\dot{\omega}=-\frac{\omega}{\tau}+\frac{\Omega}{b}+\sqrt{\frac{2q}{\tau}}\xi(t)+F_{\text{ext}}+\mathcal{F}(\omega),$
(22a)
$\dot{\Omega}=-\frac{\Omega}{\tau_{1}}-ab\omega+\sqrt{\frac{2q_{1}ab^{2}}{\tau_{1}}}\xi_{1}(t),$
(22b)

where the limit of no memory, in the stationary state, corresponds to $a=0$.

In order to have an explicit expression for the bound, we need to solve the
equation for the optimal $u$ which in this case is a scalar. Of course, we
have to do this optimization for $\mathcal{B}_{\text{max}}$ and not for
$\mathcal{B}$ so that the optimal $u$ will be a function of
$\max(|\mathcal{F}^{\prime}|)$ and not of
$\langle\mathcal{F}^{\prime}\rangle$:

$\tilde{u}=u/\langle\omega\rangle=T^{-1}z/\langle\omega\rangle=\left(\frac{\tau^{-1}+\max(|\mathcal{F}^{\prime}|)}{bB_{0}^{2}}-\frac{ab}{\tau_{1}B_{1}^{2}}\right)\left(\frac{1}{\tau_{1}^{2}B_{1}^{2}}+\frac{1}{b^{2}B_{0}^{2}}\right)^{-1}.$
(23)

We now define the non-linear forces of interest:

$\mathcal{F}_{\text{Coulomb}}(\omega)=\begin{cases}-\omega/\tau_{\mu}\quad\text{if}\quad|\omega|\leq\mu\\\
-(\mu/\tau_{\mu})\textrm{sign}(\omega)\quad\text{if}\quad|\omega|>\mu,\end{cases}$
(24)

which describes the Coulomb friction law with a viscous regime for small
velocities $|\omega|\leq\mu$. The latter is generally needed both for
numerical reasons (i.e. to avoid discontinuous forces) but it has also a
physical motivation in frictional contact models.

The PT friction force reads:

$\mathcal{F}_{\text{PT}}(\omega)=\begin{cases}-\omega/\tau_{\mu}\quad\text{if}\quad|\omega|\leq\mu\\\
-(\mu/\tau_{\mu})\textrm{sign}(\omega)\left[1+\ln(|\omega|/\mu)\right]\quad\text{if}\quad|\omega|>\mu.\end{cases}$
(25)

In both models, $\tau_{\mu}$ and $\mu$ are positive constants. Note that Eq.
(25) coincides with Eq. (1) of [19], with the constants defined in such a way
as to have a continuous function with a continuous derivative.

In order to have an explicit expression of the bound we need the maximum
modulus of the non-linear forces, which is in both cases $1/\tau_{\mu}$

$\max(|\mathcal{F}_{\text{Coulomb}}^{\prime}|)=\frac{1}{\tau_{\mu}},\quad\max(|\mathcal{F}_{\text{PT}}^{\prime}|)=\frac{1}{\tau_{\mu}}.$
(26)

We have the following bounds for the cases with memory:

$\displaystyle\lim_{t\to\infty}\frac{1}{t}\langle\Delta\theta(t)^{2}\rangle$
$\displaystyle\geq\frac{1}{\tilde{\mathcal{B}}_{\text{max}}}\geq\frac{1}{\tilde{\mathcal{B}}_{\text{max,2}}},$
(27) $\displaystyle\tilde{\mathcal{B}}_{\text{max}}$
$\displaystyle=\left(\frac{1}{\tau}+\max(|\mathcal{F}^{\prime}|)-\frac{\tilde{u}}{b}\right)^{2}B_{0}^{-2}+\left(ab+\frac{\tilde{u}}{\tau_{1}}\right)^{2}B_{1}^{-2},$
(28) $\displaystyle\tilde{\mathcal{B}}_{\text{max,2}}$
$\displaystyle=\left(\frac{1}{\tau}+\max(|\mathcal{F}^{\prime}|)\right)^{2}B_{0}^{-2}+\left(ab\right)^{2}B_{1}^{-2},$
(29)

where $\max(|\mathcal{F}^{\prime}|)$ is defined in Eq. (26) for the two
forces, and $\tilde{u}$ is given in Eq. (23). The first inequality in Eq. (27)
represents one of our main results, while the second one is considered to show
how the optimized $\tilde{u}$ can substantially improve the bound.

###  3.1 Coulomb friction with and without memory

In Fig. 1 we show the case where both variables are at temperature
$q=q_{1}=1$, i.e. the system is at equilibrium unless the external force is
switched on. Several curves in the main graph illustrate the behavior of the
mean squared displacement (rescaled by time) for different values of the
memory coupling $a$, where we remind that $a=0$ corresponds to a lack of
memory. In the inset (black dots), we consider the diffusion coefficient
$2D=\lim_{t\to\infty}\langle[\theta(t)-\theta(0)]^{2}\rangle/t$, as a function
of $a$. The first observation is the non-monotonic behaviour of the diffusion
coefficient as a function of $a$. In the descending regime (the values of $a$
for which $D$ decreases with $a$) the curve of $MSD/t$ appears to be also non-
monotonous. Apparently, the second time-scale introduced by the memory is
responsible for a kind of dynamical slowing down which eventually influences
the asymptotic diffusion, reducing $D$. In the inset, we also report the two
bounds $1/B_{\text{max}}$ and $1/B_{\text{max,2}}$. The tightest bound fairly
approximates the diffusion coefficient, including its non-monotonic behaviour.
The non-optimised bound, on the contrary, gets close to $D$ only in the limit
$a\to 0$ (no memory), while in the rest of the range it significantly
underestimates $D$, showing a monotonic trend.

Figure 1: MSD, diffusion coefficient and analytical bounds for the Coulomb
friction model with $q=q_{1}=1$. Other parameters are: $\tau=1$, $b=0.01$,
$F_{\text{ext}}=0.1$, $\tau_{\mu}=0.1$, $\mu=1$, $\tau_{1}=0.1$. Here and in
the following, we measure time in $\tau$ units. The integration time-step used
in the numerical simulations is $dt=10^{-3}$.

In Fig. 2 we repeat the same analysis in a weakly out-of-equilibrium case
$q=1.5$ and $q=1$ (the variable representing memory is slightly hotter than
the main velocity variable). The scenario looks similar to the equilibrium
case, with a non-monotonic behavior of both $MSD/t$ vs. $t$ (for large $a$)
and of $D$ vs $a$. In the inset, we report the two bounds $1/B_{\text{max}}$
and $1/B_{\text{max,2}}$, finding the same behaviour observed in the previous
case.

Figure 2: MSD, diffusion coefficient and analytical bounds for the Coulomb
friction model with $q=1.5$ and $q_{1}=1$. Other parameters as in Fig. 1.

Last, we consider a strongly out-of-equilibrium case with $q=10$ and
$q_{1}=1$, in Fig. 3. Here, while the curves of $MSD/t$ vs. $t$ are still non-
monotonous when memory is strong (large $a$), the scenario for $D$ vs. $a$ is
different from the previous cases, as it is monotonously decreasing. The bound
is confirmed and also monotonously decreasing, but it heavily underestimates
$D$ in the limit of weak memory.

Figure 3: MSD, diffusion coefficient and analytical bounds for the Coulomb
friction model with $q=10$ and $q_{1}=1$. Other parameters as in Fig. 1.

###  3.2 Prandtl-Thomlinson friction with and without Memory.

In this Section, we study the PT frictional model. Figures 4, 5 and 6 show the
results of simulations for equilibrium, weakly out-of-equilibrium and strongly
out-of-equilibrium cases, respectively. The scenario is quite similar to that
observed in the Coulomb model. At large memory the $MSD/t$ vs $t$ curve is
non-monotonous, and $D$ vs. $a$ is also non-monotonous in the equilibrium and
weakly out-of-equilibrium regimes. In these two regimes the optimised bound
works very well, and much better than the non-optimised case. In the strongly
out-of-equilibrium case, the non-monotonicy of $D$ vs. $a$ disappears and the
two bounds become indistinguishable. Remarkably, for this model, the bound
provides a nice approximation of $D$ for a much larger range of memory
coupling values $a$, with respect to the Coulomb case.

Figure 4: MSD, diffusion coefficient and analytical bounds for the PT friction
model with $q=q_{1}=1$. Other parameters are: $\tau=1$, $b=0.01$,
$F_{\text{ext}}=0.1$, $\tau_{\mu}=0.1$, $\mu=1$, $\tau_{1}=0.1$.  Figure 5:
MSD, diffusion coefficient and analytical bounds for the PT friction model
with $q=1.5$ and $q_{1}=1$. Other parameters as in Fig. 4 Figure 6: MSD,
diffusion coefficient and analytical bounds for the PT friction model with
$q=10$ and $q_{1}=1$. Other parameters as in Fig. 4

###  3.3 Discussion of analytical results for the gradient cases

The non-linear friction models considered here are, usually, not analytically
solvable. Nevertheless, in the absence of memory, the single variable
stochastic differential equation becomes a gradient system and this gives
easily an expression for the steady-state probability distribution:

$P(\omega)=\frac{\mathcal{N}\tau}{2q}\exp\left[\frac{1}{q}\int_{0}^{\omega}d\omega^{\prime}\left(-\omega^{\prime}+\tau
F_{\text{ext}}+\tau\mathcal{F}_{\text{Cou}}(\omega^{\prime})\right)\right],$
(30)

where $\mathcal{N}$ is the normalization factor. The above formula is useful
to get analytical results for the quantities $\mathcal{B}$ and $\mathcal{I}$
which are included in the most general bound expression, Eq. (14).

For the sake of simplicity, we consider only the Coulomb model without memory,
that is Eq. (22) with $a=0$ and with the force given by Eq. (24). We analyse
the effect of $\tau$ on the diffusion coefficient. As a preliminary analysis,
we verified that at small $\tau$, where the force-velocity curve is almost
linear, the steady state pdf is quasi-Gaussian, while at large $\tau$, where
the force non-linearity is important, the steady distribution is substantially
non-Gaussian (see Fig. 7).

Figure 7: Stationary probability distribution computed from Eq. (30) and from
numerical simulations of the model compared with Gaussian distributions having
the same mean and variance. Coulomb friction model without memory and $q=1$,
$\tau_{\mu}=0.5$, $\mu=1$, $F_{\text{ext}}=1.3$. For $\tau=0.1$ (panel a), we
have a quasi-Gaussian pdf, for $\tau=1.3$ (panel b) deviations from
Gaussianity are more evident.

We also take into account the effect of an external force $F_{\text{ext}}$,
which is necessary in order to compare our results with the standard TUR,
since its expression requires $F_{\text{ext}}\neq 0$. Indeed, contrary to what
happens for our bound in Eq. (14), where we have simplified the terms
$\langle\omega\rangle^{2}$, in the bound derived from the standard TUR Eq.
(18) the limit $F_{\text{ext}}\to 0$ can be singular depending on the implicit
dependence of $\langle\omega\rangle^{2}$ and $\langle\dot{S}\rangle$ on
$F_{\text{ext}}$.

Figure 8: Diffusion coefficient and analytical bounds as a function of $\tau$
for the Coulomb friction model without memory and $q=1$, $\tau_{\mu}=0.5$,
$\mu=1$.

In Fig. 8 we show the behavior of the fully analytical bound as a function of
$\tau$, for a case with a vanishing external force (panel a) and a non-zero
external force (panel b). For this model, the limit $F_{\text{ext}\to 0}$ of
$\langle\omega\rangle^{2}/\langle\dot{S}\rangle$ is regular so we can do the
comparison with the bound coming from the standard TUR in panel a. The blue
curves correspond to the diffusion coefficient measured from numerical
simulations, the black dashed lines to the bound obtained from the standard
TUR, the red solid lines correspond to our bound discussed here and the green
dash-dotted curves correspond to our bound where
$\langle\mathcal{F^{\prime}}\rangle$ is replaced with
$\max(|\mathcal{F}^{\prime}|)$. We first note that the non-monotonous behavior
in $\tau$ is well captured by all the bounds. Then, we underline that when
$F_{\text{ext}}\sim 0$ the standard TUR usually surpasses the TUR discussed
here. The opposite occurs when an external force is present, particularly the
TUR discussed here is largely dominating in the region, at large $\tau$, where
non-linear and non-Gaussian effects are important.

Figure 9: Diffusion coefficient and analytical bounds as a function of
$F_{\text{ext}}$ for the Coulomb friction model without memory and $\tau=1.3$,
$q=1$, $\tau_{\mu}=0.5$.

In Fig. 9 we show the behavior of the fully analytical bound as functions of
the external force $F_{\text{ext}}$. In panel a), the effect discussed
previously is enhanced: the TUR derived in this paper gives a tighter bound
with respect to the one obtained from the standard TUR for
$F_{\text{ext}}\gtrsim 0.6$. In panel b, we show the results for regular
Brownian diffusion. As we expect, our bound and the one coming from the
standard TUR coincide in this case.

##  4 Conclusions

We have reconsidered the information-based derivation of the Thermodynamic
Uncertainty Relations, recognising that for a class of Langevin systems with
linear coupling and a non-linear friction in the first (velocity) variable, a
new thermodynamic uncertainty relation can be obtained. The main application
of linearly coupled Langevin equations is for non-Markovian physical models
with exponentially decaying memory kernels: each auxiliary variable in the
system represents a particular memory time-scale. We remind that it is
possible to represent a power-law decaying memory kernel as a sum of a number
of exponentially decaying memory kernels (i.e. a number of auxiliary
variables), provided that such a number is sufficiently large (typically one
for each decade of the power law).

The new uncertainty relation is a bound from below for the mean squared
displacement as a function of time. It implies a bound for the diffusion
coefficient, a quantity which - in cases with non-linear forces and in
particular when more variables are present - is difficult to estimate from the
knowledge of the model parameters or from empirical data of other quantities.
The analytical bound requires the knowledge of information-theory-related
quantities. In order to simplify the expression of the bound, we also show
less tight bounds where simpler quantities are required.

We notice that, in some cases, the standard TUR gives trivial or useless
bounds when the external force is absent. The new TUR, on the contrary, due to
the fact that the entropy production related to the thermal gradients in the
system is absent from the denominator, allows one to obtain bounds (not far
from the correct value) for the diffusion coefficient also when the external
force is absent. When the external force is present we have verified that the
TUR derived in this paper can become tighter than the standard one when the
system is strongly out of equilibrium.

##  5 Acknowledgments

The authors acknowledge useful discussions with Benjamin Walter. Andrea
Puglisi and Alessandro Sarracino also acknowledge the Isaac Newton Institute
for Mathematical Sciences for support and hospitality during the programme
‘Mathematics of Movement: an interdisciplinary approach to mutual challenges
in animal ecology and cell biology’, when part of the work on this paper was
undertaken, supported by the EPSRC Grant Number EP/R014604/1.

## References

  * [1] A. C. Barato and U. Seifert.  Thermodynamic uncertainty relation for biomolecular processes.  Phys. Rev. Lett., 114:158101, Apr 2015.
  * [2] A. Baule, H. Touchette, and E. Cohen.  Stick–slip motion of solids with dry friction subject to random vibrations and an external field.  Nonlinearity, 24(2):351, 2010.
  * [3] I. Cantat.  Liquid meniscus friction on a wet plate: Bubbles, lamellae, and foams).  Physics of Fluids, 25(3):031303, 03 2013.
  * [4] L. Cerino and A. Puglisi.  Entropy production for velocity-dependent macroscopic forces: The problem of dissipation without fluctuations.  Europhysics Letters, 111(4):40012, 2015.
  * [5] A. Dechant.  Multidimensional thermodynamic uncertainty relations.  Journal of Physics A: Mathematical and Theoretical, 52(3):035001, 2018.
  * [6] C. Dieball and A. Godec.  Direct route to thermodynamic uncertainty relations and their saturation.  Physical Review Letters, 130(8):087101, 2023.
  * [7] P. G. d. Gennes.  Brownian motion with dry friction.  Journal of Statistical Physics, 119:953–962, 2005.
  * [8] T. R. Gingrich, J. M. Horowitz, N. Perunov, and J. L. England.  Dissipation bounds all steady-state current fluctuations.  Phys. Rev. Lett., 116:120601, Mar 2016.
  * [9] D. Hartich and A. c. v. Godec.  Thermodynamic uncertainty relation bounds the extent of anomalous diffusion.  Phys. Rev. Lett., 127:080601, Aug 2021.
  * [10] Y. Hasegawa and T. Van Vu.  Uncertainty relations in stochastic processes: An information inequality approach.  Phys. Rev. E, 99:062126, Jun 2019.
  * [11] H. Hayakawa.  Langevin equation with coulomb friction.  Physica D: Nonlinear Phenomena, 205(1-4):48–56, 2005.
  * [12] J. M. Horowitz and T. R. Gingrich.  Thermodynamic uncertainty relations constrain non-equilibrium fluctuations.  Nature Physics, 16(1):15–20, 2020.
  * [13] W. Hwang and C. Hyeon.  Energetic costs, precision, and transport efficiency of molecular motors.  The journal of physical chemistry letters, 9(3):513–520, 2018.
  * [14] D.-K. Kim, Y. Bae, S. Lee, and H. Jeong.  Learning entropy production via neural networks.  Physical Review Letters, 125(14):140604, 2020.
  * [15] T. Lequy and A. M. Menzel.  Stochastic motion under nonlinear friction representing shear thinning.  Physical Review E, 108(6), Dec. 2023.
  * [16] S. A. Loos.  Stochastic systems with time delay: probabilistic and thermodynamic descriptions of non-Markovian processes far from equilibrium.  Springer Nature, 2021.
  * [17] A. Manacorda, A. Puglisi, and A. Sarracino.  Coulomb friction driving brownian motors.  Communications in Theoretical Physics, 62(4):505, 2014.
  * [18] U. M. B. Marconi, A. Puglisi, L. Rondoni, and A. Vulpiani.  Fluctuation–dissipation: response theory in statistical physics.  Physics reports, 461(4-6):111–195, 2008.
  * [19] M. H. Müser.  Velocity dependence of kinetic friction in the prandtl-tomlinson model.  Phys. Rev. B, 84:125419, Sep 2011.
  * [20] A. Plati, A. Puglisi, and A. Sarracino.  Thermodynamic bounds for diffusion in nonequilibrium systems with multiple timescales.  Physical Review E, 107(4):044132, 2023.
  * [21] A. Plyukhin and A. Froese.  Nonlinear dissipation effect in brownian relaxation.  Physical Review E, 76(3):031121, 2007.
  * [22] A. Puglisi and D. Villamaina.  Irreversible effects of memory.  EPL (Europhysics Letters), 88(3):30004, 2009.
  * [23] A. Sarracino.  Time asymmetry of the kramers equation with nonlinear friction: Fluctuation-dissipation relation and ratchet effect.  Physical Review E, 88(5):052124, 2013.
  * [24] A. Sarracino, A. Gnoli, and A. Puglisi.  Ratchet effect driven by coulomb friction: The asymmetric rayleigh piston.  Physical Review E, 87(4):040101, 2013.
  * [25] U. Seifert.  Stochastic thermodynamics: From principles to the cost of precision.  Physica A: Statistical Mechanics and its Applications, 504:176–191, 2018.
  * [26] M. K. Steven.  Fundamentals of statistical signal processing.  PTR Prentice-Hall, Englewood Cliffs, NJ, 10(151045):148, 1993.
  * [27] T. Van Vu and Y. Hasegawa.  Uncertainty relations for underdamped langevin dynamics.  Phys. Rev. E, 100:032130, Sep 2019.
  * [28] A. Vanossi, N. Manini, M. Urbakh, S. Zapperi, and E. Tosatti.  Colloquium: Modeling friction: From nanoscale to mesoscale.  Reviews of Modern Physics, 85(2):529, 2013.