Rice University, United States of
AmericaHung.D.Pham@rice.eduhttps://orcid.org/0000-0002-1825-0097 Rice
University, United States of<EMAIL_ADDRESS>Indian
Institute of Technology Hyderabad<EMAIL_ADDRESS>Duc Hung Pham,
Krishna V. Palem and M.V. Panduranga Rao[100]Theory of computation Design and
analysis of algorithms [100]Theory of computationThis work was supported in
part by a Guggenheim fellowship and by DARPA under contract number
FA8750-16-2-0004.John Q. Open and Joan R. Access 2 42nd Conference on Very
Important Topics (CVIT 2016) CVIT 2016 CVIT 2016 December 24–27, 2016 Little
Whinging, United Kingdom 42 23

# A polynomial time parallel algorithm for graph isomorphism using a
quasipolynomial number of processors

Duc Hung Pham  Krishna V. Palem  M. V. Panduranga Rao

###### Abstract

Graph Isomorphism (GI) problem is a theoretically interesting problem because
it has not been proven to be in P nor to be NP complete. It is known that the
GI problem is in the low hierarchy of NP, and does not equal NP unless the
polynomial-time hierarchy collapses to some finite level. For thirty years,
the best algorithm for GI had sub-exponential running time, until Babai made a
breakthrough in 2015 when announcing a quasipolynomial time algorithm for GI
problem. Babai work gives the most theoretically efficient algorithm for GI,
as well as a strong evidence favoring the idea that class GI $\neq$ NP and
thus P $\neq$ NP. Based on Babai’s algorithm, we prove that GI can further be
solved by a parallel algorithm that runs in polynomial time using a
quasipolynomial number of processors. We achieve that result by identifying
the bottlenecks in Babai’s algorithms and parallelizing them. In particular,
we prove that color refinement can be computed in parallel logarithmic time
using a polynomial number of processors, and the $k$-dimensional WL refinement
can be computed in parallel polynomial time using a quasipolynomial number of
processors. Our work suggests that Graph Isomorphism and GI-complete problems
can be computed efficiently in a parallel computer, and provides insights on
speeding up parallel GI programs in practice.

###### keywords:

graph isomorphism, algorithm, parallel, quasipolynomial, WL refinement

###### category:

##  1 Introduction

###  1.1 Overview

An isomorphism between two graphs $G_{1}=(V_{1},E_{1})$ and
$G_{2}=(V_{2},E_{2})$ is a bijection $f:V_{1}\rightarrow V_{2}$ such that for
every $e=(u,v)\in E_{1}$, $e^{\prime}=(f(u),f(v))$ is in $E_{2}$. The Graph
Isomorphism (GI) problem is the problem that given two graphs $G_{1},G_{2}$,
determine if there exists an isomorphism between them.

The Graph Isomorphism problem is an interesting problem. Even though it is in
NP, researchers have not successfully proven it to be in P, or to be NP-
complete. In fact, it is widely believed to lie in the NP-intermediate class
(NPI), the non-empty class which lies in NP but outside P in the case that P
$\neq$ NP [8].

The Graph Isomorphism problem have witnessed advancement in the practical
aspect more than in the theoretical aspect. In practice, notable programs
include Nauty [12] and Traces [13] by McKay, saucy by Darga [4], conauto by
Presa [10], and bliss by Junttila [6]. The practical GI solvers can solve GI
of random graphs of about 10000 vertices quickly [12]. In theory, however, the
best known general theoretical bound for GI problem before Babai’s work was
sub-exponential time $2^{O(\sqrt{nlogn}})$. The algorithm was developed by
Luks and Babai in [2] in 1983. It took more than three decades for the next
best bound, quasipolynomial, to be developed. Intuitively, Babai’s work pushes
GI closer to the class of P, providing more evidence suggesting that GI is not
NP-complete. While the algorithm’s quasipolynomial runtime is still
superpolynomial, it is now remarkably closer to polynomial than to exponential
time.

On the other hand, demand for parallel algorithms is increasing. Nowadays many
computational tasks are connected with big data or large input. On the other
hand, more and more computers are equipped with multiple cores, and we have
supercomputers with millions of cores. While memory access is still a big
issue in parallel computing, it is undeniable that there is an increasingly
high demand for parallel algorithms.

###  1.2 Result

In this work, we prove that we can solve the Graph Isomorphism problem in
parallel polynomial time using a quasipolynomial number of processors. Using
Babai’s algorithm as the basis of parallelization, we demonstrate that Babai’s
algorithm, as well as GI algorithms in general, are highly parallel, and can
be sped up superpolynomially in a parallel computer.

######  Theorem 1.1.

The Graph Isomorphism problem can be solved using a parallel polynomial time
algorithm with a quasipolynomial number of processors.

Later in this paper, we will show that Babai’s work, as well as GI algorithms
in general, uses a backtracking tree structure. The tree structure enables
simple parallelization scheme. However, in Babai’s algorithm the computation
of each node in the tree is superpolynomial, and therein lies challenge of
parallelizing the tree nodes. Several operations in the computation of each
node that are bottlenecks and we will develop a parallelization scheme for
each of them. One particular nontrivial parallelization is the parallelization
of the $k$-dimensional WL refinement procedure, which is a highly sequential
procedure where each iteration depends on the result of the previous
iteration.

The structure of the paper is as follows. Section 2 presents preliminary
techniques and a short overview of Babai’s algorithm. Section 3 identifies the
bottlenecks in Babai’s algorithm and parallelization schemes designed to
handle them. Section 4 presents our conclusion and the impact of this result.

##  2 Preliminaries

###  2.1 Combinatorial techniques: Color refinement and $k$ dimensional WL
refinement, individualization

Color refinement has been used in most of the existing practical or
theoretical GI algorithms. The idea behind color refinement is to use colors
to represent similarities. Let us use $C(v)$ to denote the color of a vertex
$v$. A vertex $v$ from graph $G_{1}$ with color $C(v)$ can only be mapped to a
vertex $u$ in $G_{2}$ with color $C(u)$ if $C(u)\equiv C(v)$. Similarity can
be evaluated using multiple criteria, for example the degree of a vertex. A
common color refinement scheme is to evaluate the similarity of two vertices
based on the colors of their neighbors, i.e. $C(u)\neq C(v)$ if
$\\{C(u^{\prime})|u^{\prime}\in
neighbors(u)\\}\neq\\{C(v^{\prime})|v^{\prime}\in neighbors(v)\\}$. After
comparing the vertices and updating the new colors, the process repeats until
equilibrium is achieved - no additional color class is created. Two graphs are
isomorphic only if they have the same coloring after refinement. Note that we
call the coloring, as well as any partition scheme, canonical if they must be
preserved by isomorphism and therefore must be agreed between two graphs.

Color refinement is a powerful tool as the GI problem for random graphs can
usually be solved using one iteration of color refinement [3]. However, color
refinement fails to distinguish graphs in hard cases, a notable example being
the case of regular graphs. Therefore, an extension of color refinement has
been developed by Weisfeiler-Lehman [15]: WL refinement, and then later
$k$-dimensional WL refinement. Note that color refinement is usually referred
to as 1-dimensional WL refinement.

The classic WL refinement is an extension of color refinement. For a graph
$G=(V,E)$ and two elements $x,y\in V$, the new color $C(x,y)$ encodes the old
color and, for all $j,k$, the number of elements $z\in\Omega$ such that
$C(x,z)=j$ and $C(y,z)=k$. Thus, the classic WL refinement refines the colors
of not only the vertices but also the edges. From construction, the classic WL
refinement is also called the 2-WL refinement.

The $k$-dimensional WL refinement further leverages the classic technique to
work on tuples of size $k$ of the vertices. For a $k$-tuple $(x_{1},\ldots
x_{k})$ of vertices, the new color encodes the old color and for all tuples
$\\{j_{1},j_{2},\ldots j_{k}\\}$, the number of elements $y$ such that
$(\forall i<k)$ $C(x_{1},\ldots,x_{i-1},y,x_{i+1},\ldots,x_{k})=j_{i}$. It is
worth noting that the classic WL refinement and the $k$-dimensional WL
refinement can be extended to work with relational structures, a
generalization of graphs, as in Babai’s paper [1].

On the other hand, individualization is a technique that can complement color
refinement. Given two graphs $G_{1}$ and $G_{2}$, individualization maps a
fixed vertex $v$ from $G_{1}$ to a vertex in $G_{2}$, usually by giving them a
unique color. Since multiple vertices in $G_{2}$ could potentially be mapped
to $v$ (i.e. have the same color as $v$), a GI algorithm that uses
individualization has to consider all possible mappings. Therefore,
individualization generates a multiplicative cost to the runtime of the
program. Combining individualization and color refinement will produce a
program with a backtracking search tree structure, see McKay [12].

###  2.2 Group theoretic techniques: Group procedures

Group theory appears naturally in the course of solving Graph Isomorphism,
because the set of automorphisms of a graph $G$ naturally makes up a group,
and the set of isomorphisms between two graphs $G$ and $G^{\prime}$ must be a
coset of that group.

Here we present some basic permutation group concepts and notations. A
permutation group is a group that contains the permutations of a set of
objects. The set of objects is called the permutation domain, or domain set,
and we say that the permutation group acts on that domain. The symmetric group
of a domain $\Omega$ is the group containing all permutations of that domain,
denoting $S_{\Omega}$ or $S_{n}$ if $|\Omega|=n$. Given a permutation group
$\mathcal{G}$ subgroup of $S_{\Omega}$ (denotes $\mathcal{G}\leq S_{\Omega}$),
a $\mathcal{G}$ orbit of a set element $a\in\Omega$ is the set of elements to
which $a$ can be mapped through the group action. Orbits therefore are
equivalence classes of objects that can be mapped to each other. If a group
has only one orbit, then the group is called transitive i.e. every pair of
objects can be mapped to each other. If a group acts transitively on a domain
set, a block is a subset that either gets mapped to itself, or gets translated
somewhere else entirely. The block system naturally forms a partition of the
domain set, and the minimal block system is the block system with the minimum
number of blocks (the block size is maximum).

In GI computation, a group is usually represented using its set of generators.
We will list some useful computation group operations that appear in the
literature. Given a group $\mathcal{G}\leq S_{n}$) with generator set $A$, we
can

  * •

find the orbits of $\mathcal{G}$ and transitivity of $\mathcal{G}$ by
iterating through the images of the domain set elements.

  * •

find $\mathcal{G}$’s minimal block system by following [5], who quotes [11]
and [7]: fix $a\in\Omega$, for $b\neq a\in\Omega$, construct the graph with
$\Omega$ its set of vertices and $\\{\\{a,b\\}^{g},g\in\mathcal{G}\\}$ as
edges, then the connected component containing $a$ and $b$ is the smallest
block containing $a$ and $b$. The action of $\mathcal{G}$ is imprimitive if
and only if the constructed graph is not connected for an arbitrary $a$ and at
least one $b$, thus we obtain the block containing $\\{a,b\\}$ thus the block
system must be non-trivial.

  * •

generate the "whole" group given a set of generators by using the Furst-
Hopcroft-Luks (FHL) algorithm [11]. From FHL, we can also determine whether an
element belongs to a group, and determine the subgroup of a group given a
polynomial time membership testing of the subgroup.

It is important to note that all these operations run in polynomial time in
$n$ and $|A|$.

###  2.3 Parallel computing and the PRAM model

In this work, we build our parallel algorithm using the Parallel RAM (PRAM)
computational model. Unlike the typical Random Access Memory (RAM) model where
there is one processor and instructions are executed sequentially, the PRAM
model features multiple processors, local memory of each processor and global
shared memory, and free read and write accesses to the global shared memory.

###  2.4 Babai’s quasipolynomial time Graph Isomorphism algorithm

In 2015, Babai showed that the String Isomorphism problem, a more general
extension of GI, can be solved in quasipolynomial time. The String Isomorphism
(SI) problem, given two string inputs $x$ and $y$ as functions from indices
set $\Omega$ to symbols $\Sigma$ and a subgroup $\mathcal{G}\leq S_{\Omega}$,
determines whether or not there exists at least one element
$\tau\in\mathcal{G}$ such that $\tau$ is an isomorphism between $x$ and $y$,
i.e. $\forall i\in\Omega$ $x(i)=y(\tau(i))$. The Graph Isomorphism problem can
be reduced to the String Isomorphism problem by flattening the input graphs’
adjacency matrices into the strings and constructing $\mathcal{G}$
corresponding to this transformation.

######  Theorem 2.1 (Babai).

The String Isomorphism problem can be solved in quasipolynomial time. As a
result, the Graph Isomorphism problem can be solved in quasipolynomial time.

In the interest of brevity, we will go through the novel concepts that make
the algorithm succeed and will give a high level description of the algorithm.
We refer the readers to the work of Babai [1] and Helfgott’s explanatory
document [5] to know more about the details of the algorithm.

Babai’s algorithm includes several key features that combine and complement
each other beautifully to make the algorithm succeed:

  * •

Divide and conquer, and recursion: at almost any point in the algorithm, the
high-level goal is to either recursively reduce the underlying group into a
combination of its subgroup and cosets, or the domain set into smaller
subsets, and consider smaller problems having those smaller subgroups or
subsets as input. This continues until the input sizes become sufficiently
small such that each problem instance can be solved in polynomial time using
brute force. The algorithm runs in quasipolynomial time if at any recursive
step the problem can be broken down into a quasipolynomial number of smaller
instances, and the depth of recursion is at most polylogarithmic.

  * •

Group theoretic and combinatorial techniques: Both combinatorial techniques
and group theoretic techniques are used in Babai’s algorithm. Intuitively,
combinatorial techniques, namely partitioning techniques, work in a top-down
fashion - those techniques try to partition the graph based on high level
asymmetry. In contrast, group theoretic techniques approach the problem bottom
up and try to construct the automorphism group element by element. Generally,
combinatorial techniques work well when the partitioning structure is highly
asymmetric, whereas group theoretical technique work well when the
partitioning structure is highly symmetric.

  * •

Symmetry vs asymmetry: Both symmetry and asymmetry benefit the algorithm.
Asymmetry enables combinatorial techniques to reduce the problem to smaller
instances as stated above. On the other hand, symmetry enables canonical
splitting of the domain set into the symmetric part and the asymmetric part,
where the automorphism group of the symmetric part is easy to construct due to
its structure.

There are five main procedures in Babai’s algorithm: Luks’ reduction, Local
Certificate, Aggregate Certificates, Design Lemma and Split-or-Johnson. Luks’
reduction is the first step in each recursion, where Babai employs the
framework developed by Luks [11] and reduces the bottleneck of GI to a certain
type of group, whose structure enables the transformation of $\Omega$ to an
auxiliary domain $\Gamma$. Local Certificate examines a logarithmic-size
subset of $\Gamma$ to determine whether it is highly symmetric with respect to
the current underlying group. Aggregate Certificates combines the results of
Local Certificates of all those subsets to infer properties about the global
symmetry or asymmetry. If there are a lot of symmetric test sets, then the
global symmetries allow efficient recursion. Otherwise, if there are a lot of
asymmetric test sets, the graph must be highly asymmetric. Design Lemma takes
a highly asymmetric relational structure as the input and produces a
uniprimitive coherent configuration, which is an important algebraic structure
closely related to the 2-dimensional WL refinement, in the worst case.
Finally, Split-or-Johnson takes a uniprimitive coherent configuration and
produces either a "split" i.e. a canonical partition of the domain set, or a
large canonically embedded graph whose structure is well-known and thus
enables efficient recursion.

##  3 Parallel Graph Isomorphism algorithm

###  3.1 Identifying the bottlenecks

In this section, we identify all the superpolynomial bottlenecks in Babai’s
algorithm. We refer the reader to the technical details of Babai’s paper [1]
or to Helfgott’s detailed time complexity analysis of Babai’s algorithm [5] to
verify that the bottlenecks mentioned are indeed all superpolynomial
bottlenecks in the algorithm. We will also identify the bottlenecks that are
easy to solve, i.e. embarrassingly parallel ones.

In Babai’s algorithm, there are two main types of bottleneck: multiplicative
bottlenecks and large computation bottlenecks.

####  3.1.1 Multiplicative bottlenecks

Multiplicative bottlenecks are bottlenecks that arise from either
individualization steps or recursion steps. These steps all create change to
the structure of the program similar to branching out from a tree, and
therefore introduce a multiplicative cost to the program.

Individualization is a step that occurs in the Design Lemma procedure, the
Split-or-Johnson procedure and the Aggregate Certificates procedure. In fact,
individualization is the main action of the Design Lemma procedure. As stated
in subsection 2.1, individualization induces a multiplicative cost to the
program’s runtime.

Recursion or reduction of a problem instance to multiple smaller instances
also cause a multiplicative cost to the program’s runtime. The reduction steps
create smaller problem instances with no guarantee about runtime (we can only
reason about the maximum depth of recursion). The steps occur in the reduction
of the working group to its subgroup and the cosets, the reduction to a
Johnson group in Luks framework, and in the recursive call to a
quasipolynomial number of instances of SI on smaller inputs in the Local
Certificate procedure.

Due to their tree-like structure, the computations of the multiplicative
bottlenecks are embarrassingly parallel.

######  Lemma 3.1.

Multiplicative bottlenecks can be computed in parallel in polynomial time,
such that any multiplicative cost to the runtime of the programs becomes
multiplicative cost to the number of processors. Furthermore, the
parallelization is work-preserving.

This case is straightforward. For the recursion/reduction steps, we assign a
processor for each problem instance created. Similarly for the
individualization steps, we assign a processor for each of the vertex choice
possibilities. Since the recursion/reduction problem instances as well as the
possible cases of individualization are independent, they can be solved in
parallel independently. Therefore, their multiplicative cost to the runtime of
the program becomes multiplicative cost to the number of processors needed.

####  3.1.2 Large computation bottlenecks

Large computation bottlenecks are the superpolynomial tasks which appear
throughout Babai’s algorithm. They are: 1) the large repetitive tasks, 2) the
aggregating results and procedures on large groups tasks, and 3) the
$k$-dimensional WL refinement task.

Large repetitive tasks are computational tasks that consist of a process that
repeats itself multiple times independently. There are two such tasks in
Babai’s algorithm: one is the repeating of Local Certificate procedure for
$\binom{|\Gamma|}{k}$ test sets where $k=\log{}n$ at the start of Aggregating
Certificates, and the other is the search through possibly all tuples of $\leq
k$ elements to see if individualizing those elements yields the desired
outcome. Again, due to each iteration being independent, those tasks are
embarrassingly parallel.

######  Lemma 3.2.

Large repetitive tasks can be computed in parallel in polynomial time.
Furthermore, the parallelization is work-preserving.

Aggregating results tasks arise when the algorithm needs to combine results
from multiple smaller problem instances or subprocedures. This type of task
occurs in multiple places in Babai’s algorithm. For example, almost all
reduction steps, such as the reduction step in Luks framework or the reduction
step in the Local Certificate procedure, require combining the resulting
groups or cosets (represented by their generating sets) into one group by
merging their generating sets. The Aggregating Certificates procedure also
combines the certificates into one group in a similar fashion. This is not
problematic memory-wise, since we can avoid memory complications between
processors by using the PRAM model of computation. However, this creates large
group (in term of size of generating sets). As we mentioned above, many group
operations have runtime polynomial in the size of the generating set.
Therefore, this causes a superpolynomial bottleneck. In the next part of the
paper we develop a generating set refinement procedure to deal exclusively
with this type of bottleneck.

The last type of bottleneck is the $k$-dimensional WL refinement. This
operation deserves a separate subsection, since it is a highly sequential and
highly non-trivial obstacle to parallelization. Our approach to parallelize
this operation will be discussed in detail in subsection 3.3.3.

###  3.2 Group generators refinement

As stated previously, group operations usually run in polynomial time of the
size of the generating set. Therefore, operations on a group with a large
(non-polynomial) generating set is problematic.

In this subsection, we develop a subroutine that refines the generating set of
any permutation group.

######  Theorem 3.3.

Given a permutation group $G\leq S_{n}$ with generating set $A$, we can create
a new generating set $A^{\prime}$ of $G$ which has size $O(n\log{}n)$. The
process can be done in parallel in $O(n\log{}n)$ steps using $|A|$ number of
processors.

######  Proof 3.4.

The subroutine is based on the FHL algorithm (also called Schreier-Sims
algorithm), an algorithm described by Luks in [11]. Given $G$ some permutation
group with domain set $\Omega=\\{x_{1},x_{2},\ldots x_{n}\\}$, the FHL
algorithm works with a tower of stabilizers
$G=G_{0}>G_{1}>G_{2}>\ldots>G_{n-1}={e}$, where $G_{i}$ is the group of
permutations that keeps all elements $x_{1},x_{2}\ldots x_{i}$ fixed. The
algorithm builds a set of $C_{i}$ of representatives of $G_{i}/G_{i-1}$ and
thus $\bigcup_{i\leq j<n-1}C_{j}$ generates $G_{i}$ for all $i<n$.

We note that the FHL algorithm itself runs in polynomial time of the size of
the group’s generating set. Our subroutine is a derivation of the FHL
algorithm, but is parallel so that it can handle input groups with large
generating sets.

Algorithm 1 MembershipTest

Input: permutation $x$, set of coset representatives $C_{0},C_{1},\ldots
C_{n-2}$
Output: true/false, index $j$

  1. 1.

$j=0$

  2. 2.

while $j<n-1$ do:

    1. (a)

if $\gamma^{-1}x^{\prime}\in G_{i+1}$ for some $\gamma\in C_{i}$ then
$x^{\prime}\leftarrow\gamma^{-1}x^{\prime}$

    2. (b)

else return (false, $j$)

  3. 3.

return (true, -1)

Algorithm 2 RefineGeneratingSet

Input: $A=\\{x_{1},x_{2},\ldots x_{N}\\}$ a generating set of a group $G\leq
S_{n}$
Output: $A^{\prime}=\\{x^{\prime}_{1},x^{\prime}_{2},\ldots
x^{\prime}_{N^{\prime}}\\}$ a generating set of $G$ where
$N^{\prime}=O(n\log{}n)$

  1. 1.

$C_{i}\leftarrow e$ for $i=0$ to $n-2$

  2. 2.

$A^{\prime}\leftarrow\\{e\\}$

  3. 3.

while true:

    1. (a)

parfor all elements $x\in A$ do:

      1. i.

call MembershipTest($x$, $C_{1},\ldots C_{n-2}$)

    2. (b)

if all element return true, break

    3. (c)

$x\leftarrow$ a randomly chosen element from the ones that return false

    4. (d)

$A^{\prime}\leftarrow\\{A^{\prime}\cup x\\}$, $C_{j}=\\{C_{j}\cup x\\}$ for
$j$ = the index returned by MembershipTest($x$, $C_{1},\ldots C_{n-2}$)

  4. 4.

return $A^{\prime}$

Proof of correctness: Denote $G(A^{\prime})$ the group generated by generating
set $A^{\prime}$. The membership testing method is derived from FHL algorithm:
for each element $x$, for increasing $i$ the algorithm searches for a
representative of the coset of $x$ modulo $G_{i+1}$ in the set $C_{i}$. If
such a representative is found, then $x$ is already contained in the group
created by $A^{\prime}$. The RefineGeneratingSet procedure builds $A^{\prime}$
by gradually adding elements from $A$ to it gradually. In each iteration,
elements of $A$ are tested to see if they are a member of $G(A^{\prime})$. It
chooses one such element and updates $A^{\prime}$ and the $C_{i}$s. The
algorithm stops when there are no more elements in $A$ can be added to
$A^{\prime}$, in other word all elements in $A$ are contained $G(A^{\prime})$.

Proof of time complexity: First, it is easy to see that the membership testing
function runs in polynomial time of $n$ and $A^{\prime}$ (note that
$\bigcup_{i}C_{i}=A^{\prime}$). Thus, the parfor loop runs in polynomial time
using $|A|$ processors. Therefore, the runtime of the algorithm is polynomial
if the number of iterations in the main while loop (line 3), which is
$|A^{\prime}|$, is polynomial. In each of its iterations, the main while loop
finds one element that is not contained in $G(A^{\prime})$ and add it to
$A^{\prime}$. Doing this results in a $G(A^{\prime})_{new}$ that contains
$G(A^{\prime})_{old}$ as a proper subgroup, meaning $G(A^{\prime})_{new}$
contains $G(A^{\prime})_{old}$ and at least one coset. Therefore,
$|G(A^{\prime})_{new}|\geq 2|G(A^{\prime})_{old}|$ and thus
$|G(A^{\prime})|\geq 2^{k}$ where $k=|A^{\prime}|$. We note that
$G(A^{\prime})$ is always a subset of the symmetric set, thus
$|G(A^{\prime})|\leq n!$. Therefore, $2^{k}\leq n!$ which means
$k\leq\log{}(n!)<n\log{}n$.

Therefore, the newly created generating set $A^{\prime}$ has size polynomial
and the subroutine can be done in parallel polynomial time using $|A|$
processors.

###  3.3 Parallelization of $k$-dimensional WL refinement

We here introduce one of the main contributions of this paper: the
parallelization of the $k$-dimensional WL refinement.

In Babai’s algorithm, the $k$-dimensional WL refinement procedure appears in
the Split-or-Johnson procedure. Recall the definition of $k$-dimensional WL
refinement in subsection 2.1. Similar to color refinement, $k$-dimensional WL
refinement is a iterative procedure where the next iteration depends on the
results of previous iterations. No proven upperbound on the number of
iterations of $k$-dimensional WL refinement exists other than the trivial
$n^{k}$ one (each iteration procedures at least one new color and there are at
most $n^{k}$ colors). Therefore, the $k$-dimensional WL refinement procedure
is a highly non-trivial obstacle for parallelization.

In this subsection, we will derive a parallelization scheme for the
$k$-dimensional $WL$ refinement. We will first develop a parallelization
scheme for the 1-dimensional WL refinement by leveraging it to 2-dimensional
WL refinement, and then develop a transformation between performing
$k$-dimensional WL refinement and performing 1-dimensional WL refinement.

####  3.3.1 Color refinement and 2-dimensional WL refinement

Let us first consider 1-dimensional WL refinement i.e. color refinement. The
number of iterations of color refinement has a tight upperbound $O(n)$, as
there are graphs where color refinement requires a linear number of iterations
such as the line graph.

We will now take a different perspective by looking at a pair of vertices
distinguished in the $h^{th}$ iteration of color refinement. Consider a graph
$G=(V,E)$ with initial coloring $\mathcal{C}$ of the vertices. We can expand
the initial coloring to 2-tuples by assigning color of $\mathcal{C}(u,u)$ by
the color of $u$ and $\mathcal{C}(u,v)$ based on whether or not $(u,v)\in E$.
Then, we employ the idea of walk refinement by Lichter [9]. Given a $G=(V,E)$
with an initial edge coloring $C$, the $k$ walk refinement is a refinement
scheme that for $G$ and $C$ determines a new coloring $\mathcal{C}_{W[k]}$
defined by:

$\mathcal{C}_{W[k]}(u,v)=\\{\\{{\mathcal{C}}(u,w_{i_{1}},w_{i_{2}},\ldots,w_{i_{k-1}},v)\\},w_{i_{j}}\in
V\\}$

where ${\mathcal{C}}(v_{1},v_{2},\ldots
v_{k})=(C(v_{1},v_{2}),C(v_{2},v_{3})\ldots C(v_{k-1},v_{k}))$.

Intuitively, a $k$-walk refinement refines the color of a pair $(u,v)$ by
comparing the colors of all the possible walks of length $h$ from $u$ to $v$.
Note that a $k$-walk refinement also implicitly contains $k^{\prime}$-walks by
having $k-k^{\prime}$ first steps stationary.

######  Theorem 3.5.

Given a graph $G=(V,E)$ with two vertices $u$ and $v$ in $V$ distinguished in
the $h^{th}$ iteration of color refinement, $h>1$. Then $u,v$ can be
distinguished by a $2h$-walk refinement, i.e.
$\mathcal{C}_{W[2h]}(u,u)\neq\mathcal{C}_{W[2h]}(v,v)$.

######  Proof 3.6.

Proof by induction. The statement is trivial true for $h=1$ where the
algorithm simply returns the initial colors.

For $u,v$ a pair of vertices distinguished in the $h^{th}$ iterations, they
must satisfy these conditions:

the set of colors of the neighbors of $u$ = $\\{u_{1},u_{2}\ldots u_{j}\\}$
and the set of colors of the neighbors of $v$ = $\\{v_{1},v_{2},\ldots
v_{j}\\}$ (assuming they have the same number of neighbors) must be different
in the $h^{th}$ iteration.

In every iteration $h^{\prime}<h-1$, $u$ and $v$ must have exactly the same
neighbor color set.

Consider a pair of vertices $u^{\prime}$ and $v^{\prime}$ in the neighbors of
$u$ and $v$ respectively distinguished by the $(h-1)^{th}$ iteration of color
refinement. By induction hypothesis, $u^{\prime}$ and $v^{\prime}$ are also
distinguished by the $2(h-2)$ walk refinement, i.e.
$\mathcal{C}_{W[2h-2]}(u^{\prime},u^{\prime})\neq\mathcal{C}_{W[2h-2]}(v^{\prime},v^{\prime})$.

Now, consider $\mathcal{C}_{W[2h]}(u,u)$ and $\mathcal{C}_{W[2h]}(v,v)$. These
colors are distinguished using $\\{{\mathcal{C}}(u,w_{1}$ $,\ldots
w_{2h-1},u\\}$ and $\\{{\mathcal{C}}(v,w_{1},\ldots w_{2h-1},v\\}$ . We will
only consider a subset of all the walks which is $\\{u,u_{i},w_{1},\ldots
w_{2k-3},u_{i},u\\}$ and $\\{v,v_{i},w_{1},\ldots w_{2k-3},v_{i},v\\}$ for
$i=1,2\ldots j$ and $u_{1},u_{2},$ $\ldots u_{j}$ are the neighbors of $u$,
and similar for $v$. We will prove that these two subsets of walks result in a
different set of colors, and that these two subsets having different colors
result in $\mathcal{C}_{W[2h]}(u,u)\neq\mathcal{C}_{W[2h]}(v,v)$.

First, since $u$ and $v$ are distinguished in the $h^{th}$ iteration of color
refinement, their set of neighbor colors must be different after the
$(h-1)^{th}$ color refinement iteration. From the induction hypothesis, this
means that
$\\{\\{\mathcal{C}_{W[2h-2]}(u_{i},u_{i})\\}\\}\neq\\{\\{\mathcal{C}_{W[2h-2]}(v_{i},v_{i})\\}\\}$
where the $u_{i},v_{i}$ are neighbors of $u,v$ respectively. This means that
$\\{\\{{\mathcal{C}}(u_{i},w_{1},\ldots,w_{2h-3},u_{i})\\}\\}\neq\\{\\{{\mathcal{C}}(v_{i},w_{1},\ldots
w_{2h-3},v_{i})$. Therefore, ${\mathcal{C}}\\{(u,u_{i},w_{1},\ldots
w_{2k-3},u_{i},u)\\}$ $\neq$ ${\mathcal{C}}\\{(v,v_{i}$ $,w_{1},\ldots
w_{2k-3},$ $v_{i},v)\\}$, thus those two subsets of walks result in a
different set of colors.

Second, we will prove that these subsets of walks have a unique role in all
the $2h$ walks considered. This is indeed true, because $\mathcal{C}(u,u_{i})$
the initial color of an edge is unique from the definition.

Thus, if $u$ and $v$ are distinguished in the $h^{th}$ iteration of color
refinement, then they can be distinguished by a $2h$-walk refinement.

However, Lichter [9] proved that a $2h$-walk refinement can be simulated by
$\log{}(2h)$ iterations of 2-dimensional WL refinement.

######  Theorem 3.7 (Lichter).

$k$-walk refinement can be simulated with $\left\lceil{\log{}k}\right\rceil$
iterations of Weisfeiler-Leman refinements.

From the above two theorems, we immediately come to our conclusion.

######  Corollary 3.8.

Given a graph $G=(V,E)$, if color refinement stabilizes in $h$ iterations,
then we can use 2-dimensional WL refinement to simulate the result using
$\log{}2h$ iterations.

We note that this result has significant implications beyond Babai’s
algorithm.

######  Theorem 3.9.

There is a parallel algorithm to simulate the result of color refinement in
logarithmic time.

######  Proof 3.10.

We can use 2-dimensional WL refinement to simulate color refinement as stated
above. Note that each step of the 2-dimensional WL refinement is completely
parallelizable (i.e. can be done in constant time by a work-preserving
parallel algorithm), therefore the runtime of the parallel algorithm is the
number of iterations, which is $O(\log{}n)$ as from the above two theorems.

Achieving the result of color refinement in logarithmic time is remarkably
useful for practical GI solvers, which mainly use color refinement instead of
the complicated $k$-dimensional version. Furthermore, this result implies that
any decision problem which can be NC-reduced to color refinement is in NC, the
class of problems that can be solved efficiently in a parallel computer.

####  3.3.2 $k$-dimensional WL refinement and color refinement

Now we will present a way to simulate $k$-dimensional WL refinement using
color refinement

######  Theorem 3.11.

Given a graph $G=(V,E)$, there is a graph $G^{\prime}$ such that performing
$k$-dimensional WL refinement on $G^{\prime}$ can be simulated by performing
color refinement on $G^{\prime}$.

######  Proof 3.12.

Recall that the $k$-dimensional WL refinement determines the new color of a
tuple $(w_{1},w_{2}\ldots w_{k})$ based on the number of $\\{y\in V|(\forall
i<k)(C(w_{1},\ldots w_{i-1},y,w_{i+1},\ldots,w_{k})=c_{i})\\}$ for all tuples
of color $c_{1},c_{2}\ldots c_{k}$.

Denote $V$ $=\\{w_{1},w_{2},\ldots w_{n}\\}$. We construct $G^{\prime}$ as
follows: For each $k$ tuple of vertices $(w_{i_{1}},w_{i_{2}},\ldots
w_{i_{k}})$ in $G$, create a corresponding vertex $x$ in $G^{\prime}$. Let us
write $x\approx(w_{i_{1}},w_{i_{2}}\ldots w_{i_{k}})$, and call all the $x$s
the base layer. For each vertex $x$ of $G^{\prime}$ that corresponds to the
tuple $(w_{i_{1}},w_{i_{2}}\ldots w_{i_{k}})$, construct two layers of
auxiliary nodes. The first layer consists of $n$ nodes $u_{1},u_{2},\ldots
u_{n}$. Connect each of the $u_{i}$ with $x$. The second layer consists of
$nk$ nodes $v_{1,1},v_{1,2},\ldots v_{1,k},v_{2,1},v_{2,2}\ldots v_{n,k}$.
Connect each $v_{i,j}$ with $u_{i}$ and with the node that corresponds to the
$k$-tuple $(w_{i_{1}},w_{i_{2}},\ldots,w_{i_{j-1}},w_{i},w_{i_{j+1}},\ldots
w_{i_{k}})$. We say that these nodes in the auxiliary $u$ and $v$ layer come
out from $x$ in order to distinguish them with the auxiliary nodes that expand
from other vertices, or come in, to $x$, and call the layer of $u$ nodes that
expand from $x$ as constructed above the out-layer and the layer of $y$ nodes
that connect to $x$ the in-layer. The initial color of each $x$ in
$G^{\prime}$ is similar to the initial color of the corresponding $k$-tuple in
$G$, all the auxiliary $u_{i}$ are given a same new color, and for each $j$
all the auxiliary $v_{i,j}$ are given a same new color. The construction
scheme of the out-layers of node $x$ in $G^{\prime}$ corresponding to a
$k$-tuple $(w_{i_{1}},w_{i_{2}}\ldots w_{i_{k}})$ is as illustrated in the
following figure.

Figure 1: Illustration of constructing the auxiliary out-layers from a vertex
$x$ in $G^{\prime}$ to simulate performing $k$-dimensional WL refinement on
$G$

For vertex $x\in G^{\prime}$ corresponding to the tuple
$(w_{i_{1}},w_{i_{2}}\ldots w_{i_{k}})$ of vertices in $G$, denote
$\overrightarrow{W}_{j}^{y}$ as the tuple $(w_{i_{1}},\ldots
w_{i_{j-1}},y,w_{i_{j+1}},\ldots w_{k})$ and $\overrightarrow{x}_{j}^{y}$ the
corresponding vertex. Intuitively, for each $x\approx(w_{1},w_{2}\ldots
w_{k})$ the auxiliary $u_{i}$ is used to capture the set of tuples
$\overrightarrow{W}_{j}^{w_{i}}$ for all $j\leq n$, and the auxiliary
$v_{i,j}$ with distinct initial colors for different $j$ are used to make sure
that the exact ordering $j$ of the sequence of colors of
$\overrightarrow{W}_{j}^{w_{i}}$ affects the colors of $u_{i}$. Note that the
different layers are always distinguished since their initial colors are
different.

We will prove that color refinement on $G^{\prime}$ results in the same
coloring of the base layer nodes as performing $k$-dimensional WL refinement
on $G$.

Consider 2 $k$-tuples $W=(w_{i_{1}},w_{i_{2}}\ldots w_{i_{k}})$ and
$W^{\prime}=(w_{i^{\prime}_{1}},w_{i^{\prime}_{2}}\ldots w_{i^{\prime}_{k}})$
of vertices in $G$ that are not distinguished at the end of $k$-dimensional WL
refinement. In this case, the tuples’ corresponding vertices $x$ and
$x^{\prime}$ in $G^{\prime}$ must have the same out-layer. Indeed, each $i$
the $u_{i}$ nodes captures the color determined by the set of colors of
$\overrightarrow{x}_{j}^{w_{i}}$. Since $W$ and $W^{\prime}$ are not
distinguished by the $k$-dimensional WL refinement, the set of
$\overrightarrow{x}_{j}^{w_{i}}$ and $\overrightarrow{x^{\prime}}_{j}^{w_{i}}$
must be symmetric. Therefore, the $u$ nodes of $x$ and $x^{\prime}$ must be
identical up to permutation. On the other hand, $x$ and $x^{\prime}$ must also
have the same in-layer, because the color of each of the $v$ nodes coming from
$\overrightarrow{x}_{j}^{y}$ depends on the color of $x$ and of
$\overrightarrow{x}_{j}^{y}$, which must be symmetric between $x$ and
$x^{\prime}$.

Now, if $W$ and $W^{\prime}$ are distinguished at the end of $k$-dimensional
WL refinement, then it is straightforward to see that $x$ and $x^{\prime}$ are
distinguished by the end of color refinement on $G^{\prime}$. This is because
the out-layer $v_{i,j}$ of each $u_{i}$ node $x$ captures the colors of
$\overrightarrow{x}_{j}^{w_{i}}$ not preserved under permutation due to the
initial coloring of $v_{i,j}$. Therefore, asymmetries in the colors of
$\overrightarrow{x}_{j}^{w_{i}}$ and $\overrightarrow{x^{\prime}}_{j}^{w_{i}}$
(which are detected at some point because $W$ and $W^{\prime}$ are
distinguished) lead to asymmetries in the colors of the auxiliary $u_{i}$ of
$x$ and $x^{\prime}$, which in turn leads to different colors of $x$ and
$x^{\prime}$.

####  3.3.3 Parallelizing $k$-dimensional WL refinement

We arrive at our main goal of this subsection.

######  Theorem 3.13.

Given a graph $G=(V,E)$, we can perform $k$-dimensional WL refinement, where
$k=O(\log{}n)$, in parallel in polynomial time using a quasipolynomial number
of processors.

######  Proof 3.14.

First, we use the simulation in theorem 3.11 to convert performing
$k$-dimensional WL refinement on $G$ to performing color refinement on
$G^{\prime}$. From construction, $G^{\prime}$ has $O(nkn^{k})$ vertices,
therefore color refinement on $G^{\prime}$ runs at most $O(nkn^{k})$
iterations.

Next, we use 2-dimensional WL refinement on $G^{\prime}$ to simulate color
refinement on $G^{\prime}$ as in Corollary 3.8. The result of color refinement
can be simulated in $\log{(2nkn^{k})}$, which is $O(n)$ iterations of
2-dimensional WL refinement for $k=O(\log{}n)$.

Finally, we can execute this 2-dimensional WL refinement in parallel
polynomial time. We will prove that in each iteration of 2-dimensional WL
refinement, the computation can be done efficiently in parallel. Recall that
in the 2-dimensional WL refinement, the new color of a pair of vertices $u,v$
is determined by the number of elements $z\in\Omega$ such that $c(x,z)=j$ and
$c(y,z)=k$ for all $j,k$. For each pair of vertices $(u,v)$, looking at the
colors $c(u,j)$ and $c(j,v)$ for all vertices $j$ in $G^{\prime}$ is entirely
parallelizable and can be done in constant time. Comparing those sets of color
between two pair of vertices takes constant time, and performing all the
comparisons among all sets of two pairs of vertices can be done in parallel.
Therefore, each iteration of 2-dimensional WL refinement of $G^{\prime}$ can
be done in parallel polynomial time, thus whether the time complexity of a
parallel execution of 2-dimensional WL refinement is polynomial or not depends
on the total number of iterations needed, which is polynomial in $n$.
Therefore, the time complexity of simulating color refinement on $G^{\prime}$
using 2-dimensional WL refinement in parallel is polynomial, with a
quasipolynomial total amount of work (because our parallelization of
2-dimensional WL refinement iterations is work-preserving, and there is a
quasipolynomial number of vertices in $G^{\prime}$).

In conclusion, we can perform $k$-dimensional WL refinement where
$k=O(\log{}n)$ in parallel polynomial time with a quasipolynomial amount of
total work, and therefore a quasipolynomial number of processors needed.

###  3.4 The combined algorithm

We have examined all the bottlenecks in Babai’s algorithm and parallelized
them. In this subsection, we combine the parallelization schemes from the
previous parts to create the master algorithm.

######  Theorem 3.15.

The graph isomorphism problem can be solved by a parallel algorithm in
polynomial time using a quasipolynomial number of processors.

######  Proof 3.16.

From lemma 3.1 and 3.2, the multiplicative bottleneck and the large repetitive
tasks bottleneck of Babai’s algorithm can be done in parallel polynomial time
in a work-preserving manner and such that any multiplicative cost to the
runtime of the program becomes multiplicative cost to the number of
processors.

For group operations on large groups, which are formed by aggregating results
tasks, we apply our RefineGeneratingSet subroutine developed in subsection 3.2
after every aggregation of results operation. By doing this, we guarantee that
all the groups in the algorithm have a polynomial size generating set, and
therefore the runtime of group operations on them is not superpolynomial. The
parallel subroutine runs in polynomial time and introduces an additional worst
case quasipolynomial $n^{\log{}n}$ cost to the total amount of work.

For $k$-dimensional WL refinement, we follow the procedure described in
theorem 3.11 to parallelize the procedure to parallel polynomial time. The
parallelization introduces a $n^{O(\log{}n)}$ multiplicative cost to the total
amount of work. However, we note that this is an one-time cost, and thus the
total amount of work done by the algorithm after parallelization is still
quasipolynomial.

Combining the above arguments, we arrive at our parallel version of Babai’s
algorithm. The algorithm runs in polynomial time in each iteration, any
multiplicative cost to the runtime of the algorithm after each iteration is
converted into multiplicative cost to the number of parallel processors, and
there are a polynomial (in fact logarithmic) number of iterations in total.
Thus, the whole parallel algorithm runs in polynomial time. The amount of
total work is increased by at most a quasipolynomial multiplicative cost, and
thus is still quasipolynomial. Therefore, the number of processors needed is
quasipolynomial.

##  4 Conclusion and remarks

In our work, we have proven that the GI problem can be solved using a parallel
algorithm that runs in polynomial time using a quasipolynomial number of
processors. Therefore, problems in the GI complexity class, which are problems
that can be polynomially reduced to Graph Isomorphism, can be solved in
parallel polynomial time using a quasipolynomial number of processors.

This result implies that in theory solving the worst cases of the GI problem
is tractable in a parallel computer. It also implies that Babai’s algorithm is
highly parallel and can be sped up superpolynomially in a parallel computer.
We note that despite being quite complicated, Babai’s algorithm shares
multiple similarities with other GI algorithms. Many of those similarities are
bottlenecks of parallelization, such as the backtracking search tree
structure, color refinement and $k$-dimensional operations, individualization
and so on. Therefore, the parallelization techniques used in this paper can be
generalized to other GI algorithms as well, indicating that algorithms for GI
in general is highly parallel.

In addition, our finding that color refinement can be calculated in
logarithmic time in parallel can potentially have a notable impact in
practical GI solvers. On the one hand, almost all practical GI solvers use
color refinement as the main tool to solve the problem. On the other hand,
given the increase in demand of big data processing and parallel algorithms,
there have been a surprisingly few practical parallel algorithms for the GI
problem - it is hard to find one such program beside [14]. Therefore, our
techniques used for parallelizing color refinement can potentially be applied
to create new practical parallel GI algorithms, or to improve existing ones.

## References

  * [1] László Babai.  Graph isomorphism in quasipolynomial time [extended abstract].  In Proceedings of the Forty-eighth Annual ACM Symposium on Theory of Computing, STOC ’16, pages 684–697, New York, NY, USA, 2016. ACM.  URL: http://doi.acm.org/10.1145/2897518.2897542, doi:10.1145/2897518.2897542.
  * [2] László Babai and Eugene M. Luks.  Canonical labeling of graphs.  In Proceedings of the Fifteenth Annual ACM Symposium on Theory of Computing, STOC ’83, pages 171–183, New York, NY, USA, 1983. ACM.  URL: http://doi.acm.org/10.1145/800061.808746, doi:10.1145/800061.808746.
  * [3] László Babai, Paul Erdös, and Stanley Selkow.  Random graph isomorphism.  SIAM J. Comput., 9:628–635, 08 1980.  doi:10.1137/0209047.
  * [4] P. T. Darga, K. A. Sakallah, and I. L. Markov.  Faster symmetry discovery using sparsity of symmetries.  In 2008 45th ACM/IEEE Design Automation Conference, pages 149–154, June 2008.  doi:10.1145/1391469.1391509.
  * [5] Harald Andrés Helfgott, Jitendra Bajpai, and Daniele Dona.  Graph isomorphisms in quasi-polynomial time, 2017.  arXiv:arXiv:1710.04574.
  * [6] Tommi Junttila and Petteri Kaski.  Conflict propagation and component recursion for canonical labeling.  pages 151–162, 01 2011.  doi:10.1007/978-3-642-19754-3_16.
  * [7] L.G. Kovacs.  Groups - Canberra 1989: Australian National University Group Theory.  Lecture Notes in Mathematics. Springer Berlin Heidelberg, 2006.  URL: https://books.google.com/books?id=Dm96CwAAQBAJ.
  * [8] Richard E. Ladner.  On the structure of polynomial time reducibility.  J. ACM, 22:155–171, 1975.
  * [9] Moritz Lichter, Ilia Ponomarenko, and Pascal Schweitzer.  Walk refinement, walk logic, and the iteration number of the weisfeiler-leman algorithmhttps://arxiv.org/submit/2681062/preview, 05 2019.
  * [10] José López-Presa, Antonio Fernández Anta, and Luis Chiroque.  Conauto-2.0: Fast isomorphism testing and automorphism group computation.  08 2011.
  * [11] Eugene M. Luks.  Isomorphism of graphs of bounded valence can be tested in polynomial time.  In Proceedings of the 21st Annual Symposium on Foundations of Computer Science, SFCS ’80, pages 42–49, Washington, DC, USA, 1980. IEEE Computer Society.  URL: http://dx.doi.org/10.1109/SFCS.1980.24, doi:10.1109/SFCS.1980.24.
  * [12] Brendan D. McKay.  Practical graph isomorphism, 1981.
  * [13] Brendan D. McKay and Adolfo Piperno.  Practical graph isomorphism, ii.  Journal of Symbolic Computation, 60:94 – 112, 2014.  URL: http://www.sciencedirect.com/science/article/pii/S0747717113001193, doi:https://doi.org/10.1016/j.jsc.2013.09.003.
  * [14] M.-Y Son, Y.-H Kim, and B.-W Oh.  An efficient parallel algorithm for graph isomorphism on gpu using cuda.  International Journal of Engineering and Technology, 7:1840–1848, 11 2015.
  * [15] B. Yu. Weisfeiler and A. A. Leman.  Reduction of a graph to a canonical form and an algebra arising during this reduction.  1968\.


# pH modulates friction memory effects in protein folding

Benjamin A. Dalton  Freie Universität Berlin, Fachbereich Physik, 14195
Berlin, Germany  Roland R. Netz Corresponding author Freie Universität Berlin,
Fachbereich Physik, 14195 Berlin, Germany

###### Abstract

We study the non-Markovian folding dynamics of the $\alpha$3D protein under
low- and neutral-pH conditions. Recently published all-atom simulations of
$\alpha$3D by the Shaw group reveal that lowering the pH significantly reduces
both native and non-native salt-bridge interactions, which dominate the
folding dynamics. Here, we demonstrate that this physiochemical modulation
directly perturbs the folding friction, which we evaluate using non-Markovian
memory-kernel-extraction techniques. In doing so, we find that the reduction
in pH not only decreases the magnitude of the time-dependent friction acting
on the protein but also more dramatically shortens the time scale of the
friction memory effects. As a result, the folding dynamics in the low pH
system are well described by a purely Markovian model. In the neutral pH
system, however, the memory time scale is of the same order as the folding
time and is accelerated by a factor of 6 compared to a Markovian model
prediction. We demonstrate that this memory-induced barrier-crossing speed-up
is predicted by non-Markovian reaction-kinetic theories, confirming that non-
Markovian models are, in general, necessary for a quantitative description of
protein folding dynamics.

Protein folding is an essential requirement for many biological functions and
has long been the interest of theoretical and computation studies [1, 2, 3, 4,
5, 6]. To describe the folding dynamics of a protein, one typically represents
an instantaneous protein structure with a suitably chosen one-dimensional
reaction coordinate. Hence, the time evolution of the reaction coordinate is
modelled as a diffusive process over the extracted free energy profile. This
coarse-graining naturally incurs friction, representing the dissipative
interactions between the reaction coordinate and its environment. Evaluating
the friction acting on an arbitrary reaction coordinate is a non-trivial task
and depends on the model used to describe the reaction coordinate dynamics.
Markovian models have been successfully applied to model protein folding [7,
8, 9]. However, these models neglect non-Markovian effects, suggesting that
such effects are fast compared to the relevant folding time scales. Two
problems with this approach are that i) it has been shown that non-Markovian
effects can still influence barrier-crossing kinetics when memory time scales
are short compared to barrier-crossing time scales [10], and ii) proteins do
exhibit long memory time scales [11]. While non-Markovian effects have been
considered in the context of protein folding [12, 13], recent memory kernel
extraction techniques [14, 15, 16, 17] have enabled the accurate evaluation of
time-dependent friction kernels acting on arbitrary reaction coordinates and
have been validated for a range of molecular systems [18, 19, 20], thus
enabling a direct method to evaluate the relevance of non-Markovian effects.

Recently, we applied non-Markovian memory kernel extraction techniques to
study the folding and unfolding dynamics of a set of all-atom, fast-folding
protein simulations. We showed that all proteins in the study exhibit
significant non-Markovian dynamics when described using a well-characterized
fraction of native contacts reaction coordinate [11]. We showed that memory
decay times are comparable to the folding time scales and that, overall, non-
Markovian processes are not fast compared to the relevant protein folding
kinetics. Furthermore, we also showed that reaction rate theories that
explicitly incorporate non-Markovian effects are successful at predicting
folding and unfolding times across the full population of proteins. In the
present article, we apply non-Markovian techniques to study of the folding of
a small, fast-folding protein under neutral and perturbed pH conditions, and
we show that the observed deviations away from a Markovian prediction for the
folding kinetics of the neutral pH protein are as expected according to
various non-Markovian predictions for the folding kinetics of proteins.

$\alpha$3D is a designed, 73-residue, fast-folding protein [21, 22], which
folds into a tertiary structure comprised of three interacting
$\alpha$-helical subunits (Fig. 1A), and which exhibits distinct folded and
unfolded states at neutral pH, both in experiment [23] and in simulation [24].
It has been shown that the folding kinetics of the $\alpha$3D protein are
strongly influenced by non-native interactions [25, 26], in particular, by
several non-native salt-bridge interactions between positively charged lysine
and arginine side-groups and negatively charged aspartate and glutamate side
groups. These salt bridges interact between the residue-chain segments
associated with the distinct helices, even in the unfolded state [26]. The
folding kinetics of the $\alpha$3D protein were compared under low and neural
pH conditions, using both experiments and all-atom simulations [25]. It was
shown that lowering the pH significantly reduces both native and non-native
salt-bridge interactions, therefore reducing the folded state stability but,
surprisingly, increasing the folding rate. In the present article, we analyse
the all-atom simulation data for the $\alpha$3D protein in both low and neural
pH conditions (see Supplementary Materials Sec. I) and show that the pH
reduction perturbs the non-Markovian folding friction. Using the memory kernel
extraction techniques we show that both the total friction and the non-
Markovian memory time scale are significantly reduced when the pH is lowered,
so much so that the non-Markovian features that dominate the $\alpha$3D
folding kinetics at neutral pH are eliminated, allowing for a Markovian model
description for the low pH folding kinetics. The same Markovian model fails to
describe $\alpha$3D folding in neutral pH conditions, under-predicting the
folding kinetics by a factor of 6, suggesting strong non-Markovian
acceleration effects for the neutral pH system.

Figure 1:  Folding dynamics of the $\alpha$3D protein. A) Simulation snapshots
of representative states for neutral pH. All-atom trajectories taken from [24,
25]. B) 25 $\mu$s trajectory segments for the fraction of native contacts $Q$
reaction coordinate for $\alpha$3D for low pH (pH=2) and neutral pH (pH=7).
The unfolded, barrier top, and folded states are indicated by $Q_{\text{u}}$,
$Q_{\text{b}}$, and $Q_{\text{f}}$, respectively (red dashed lines). C)
Fraction of native contacts free-energy profiles $U(Q)=-k_{\rm{B}}T\log[P(Q)]$
for the two pH conditions. $P(Q)$ is the $Q$ probability density. Inset:
conditional probabilities to be on a transition path $P(TP|Q)$ at position
$Q$, suggesting that $Q$ is a good reaction coordinate for the all-atom
$\alpha$3D trajectories under both pH conditions [27]. For comparison, $Q$ is
shifted by the peak value $Q_{\text{p}}$, often identified as the transition
state, and rescaled by the distance between free energy minima
$L_{\text{U}}=Q_{\text{f}}-Q_{\text{u}}$.

To track the folding progress of the $\alpha$3D protein, we project the atomic
coordinates onto the fraction of native contacts reaction coordinate $Q$ (Fig.
1B) [28, 26] (Supplementary Materials Sec. II). While the techniques employed
in this article are, by definition, suitable for any arbitrary reaction
coordinate, irrespective of the goodness or poorness of that reaction
coordinate, we first note that $Q$ is equally good in both pH conditions,
judged by the following criteria: for the neutral pH system, $Q$ exhibits a
pronounced two-state distribution (Fig. 1C) and a sharply peaked conditional
transition-path probability $P(TP|Q)=P(Q|TP)P(TP)/P(Q)$ (Fig. 1C inset and
Supplementary Materials Sec. II), i.e. the probability to be on a transition
path when at $Q$ [27]. $P(Q)$ and $P(TP)$ are the probability density for $Q$
and the probability of being on a transition path, respectively. $P(Q|TP)$ is
the probability density conditional to being on a transition path. $P(TP|Q)$
provides a measure for the quality of a reaction coordinate, and it has been
shown that $Q$ optimally reveals the transition state for proteins [26], which
is apparent since $P(TP|Q)$ is sharply peaked and almost reaches a value of
0.5, which is an indication of a transition state in 1D Markovian theory.
These results suggest that non-Markovian effects should be negligible [13],
which we later show not to be the case.

Simulated folding times, as represented by the average residency time (the
average time spent in a state before making a transition) in the unfolded
states, of 10.2 $\pm\;2$ $\mu$s and 24.2 $\pm\;9$ $\mu$s in the low and
neutral pH conditions, respectively, can be compared to 10 $\mu$s and 20
$\mu$s, as reported previously for the same data [25], and experimental values
of $\sim 6$ $\mu$s [22] and $\sim 10$ $\mu$s [23]. Due to the strong asymmetry
in the free energy profiles (Fig. 1C) and the large differences between the
folded-state energies, we here only consider transitions which connect the
unfolded state minimum $Q_{\text{u}}$ to the barrier top $Q_{\text{b}}$.
Therefore, we define the folding times $\tau^{\text{MD}}_{\text{MFP}}$ in the
MD simulations as the mean first-passage times from $Q_{\text{u}}$ to
$Q_{\text{b}}$. Measured in this way, we evaluate folding times of
$\tau^{\text{MD}}_{\text{MFP}}=5.2\;\mu\rm{s}$ and
$\tau^{\text{MD}}_{\text{MFP}}=6.4\;\mu\rm{s}$ for the low and neutral pH
systems. These are significantly different to the residency times due to
recrossing events, i.e., excursions that reach $Q_{\text{b}}$, but which do
not continue to $Q_{\text{f}}$. In Supplementary Material Sec. III, we discuss
the various definitions for folding times and recrossing, and we present
folding time distributions and error estimates. By any definition of the
folding time, we see that the low pH system folds slightly faster than the
neutral system. As pointed out previously [25], this is interesting since the
folding barrier is slightly higher for the low pH system (Fig. 1C), so the
folding time reduction for low pH must be due to friction effects.

Figure 2: Friction and memory effects in $\alpha$3D dynamics. A) Memory
kernels $\Gamma(t)$, extracted for two different pH conditions. Inset:
magnification of the long-time behaviour. B) Running integrals $G(t)$ for the
two extracted memory kernels in A) (bold lines). The dashed lines show running
integrals constructed from four-component exponential fits to the memory
kernel data in A). Total friction $\gamma$ is given by the plateau value of
$G(t\rightarrow\infty)$, taken from the fitted curves.

To evaluate the friction, we map the $Q(t)$ trajectories for $\alpha$3D onto a
generalised Langevin equation [29, 30], given by

$\begin{split}m\ddot{Q}(t)&=-\int\limits_{0}^{t}\Gamma(t-t^{\prime})\dot{Q}(t^{\prime})dt^{\prime}-\nabla
U\big{[}Q(t)\big{]}+\xi(t).\end{split}$ (1)

$\Gamma(t)$ is the friction memory kernel, $U(Q)$ is the free energy profile
(Fig. 1C), and $m$ is the effective mass of $Q$, which is assumed to have no
position dependence [11]. $\xi(t)$ is the random force term satisfying the
fluctuation-dissipation theorem
$\langle\xi(t)\xi(t^{\prime})\rangle=k_{\rm{B}}T\Gamma(t-t^{\prime})$, where
$k_{\rm{B}}T$ is the thermal energy. Using recent memory kernel extraction
techniques, we extract the running integral of the memory kernel
$G(t)=\int_{0}^{t}\Gamma(t^{\prime})dt^{\prime}$ directly from the time series
of $Q(t)$ [16, 17], where the total friction $\gamma$ acting on $Q(t)$ is
determined by the plateau value $\gamma=G(t\rightarrow\infty)$. In Fig. 2A, we
show the memory kernels for the $\alpha$3D protein under the two pH conditions
(for additional information, see Supplementary Material Sec. IV)). Lowering
the pH in the MD simulations decreases both the amplitude and the decay time
of the memory. The latter is evident in the long-time tails, as magnified in
the inset. The memory kernel for neutral pH spans microseconds. However, for
the low pH, the decay is over an order of magnitude faster. From the running
integrals $G(t)$ in Fig. 2B, we see that the total friction $\gamma$ is
reduced by a factor of $\sim$23\. Thus, whereas the amplitude of $\Gamma(t)$
is reduced by a factor of $\sim$2 at short times, the dramatic reduction in
total friction is caused primarily by the reduced memory time.

Figure 3: Comparison of dynamic time scales for $\alpha$3D folding under two
pH conditions. A) The ratio of the memory time $\tau_{\text{mem}}$ and
diffusion time $\tau_{\text{D}}$ suggests that non-Markovian effects are more
pronounced in the neutral pH condition. B) MD simulation mean first-passage
times for transitions from unfolded state minimum to barrier top .C)
Comparison of first mean-passage times $\tau^{\text{MD}}_{\text{MFP}}$ between
$Q_{\text{u}}$ and $Q_{\text{b}}$ measured in simulations and a purely
Markovian prediction of passage times $\tau^{\text{Mar}}_{\text{MFP}}$ (Eq.
2). For the low pH condition, $\tau^{\text{MD}}_{\text{MFP}}$ is well
predicted by a Markovian model. For the neutral pH system, the transition time
is significantly faster than that measured in the MD simulation, suggesting
non-Markovian acceleration effects.

We evaluate the memory decay via the first moment of the memory kernel
$\tau_{{\text{mem}}}=\int_{0}^{\infty}t\Gamma(t)dt/\int_{0}^{\infty}\Gamma(t)dt$,
and obtain $\tau_{{\text{mem}}}=3.4\;\mu\text{s}$ and $0.07\;\mu\text{s}$ for
the neutral and low pH systems, respectively. Having extracted the friction,
we can evaluate other characteristic time scales associated with the folding
dynamics. The diffusion time
$\tau_{\text{D}}=\gamma(Q_{\text{b}}-Q_{\text{u}})^{2}/k_{\rm{B}}T$ is the
time to diffuse from the unfolded state to the barrier top in the absence of a
free energy profile (a table of all relevant time scales and extracted system
parameters is given in the Supplementary Materials Sec. I). We expect memory
effects to influence barrier crossing kinetics when
$1{\times}10^{-2}<\tau_{\text{mem}}/\tau_{\text{D}}<1{\times}10^{1}$ [10].
This condition is well satisfied for the neutral pH system, but only just so
for the low pH system (Fig. 3A). Therefore, we expect non-Markovian effects to
be more pronounced for the neutral pH system.

In Fig. 3B, we show the simulated folding mean-first passage times for the two
pH systems. To quantify deviations from Markovianity, we predict transition
times using the following Markovian overdamped formula which explicitly
accounts for the free energy profiles given in Fig. 1C [31]:

$\begin{split}\tau_{\rm{MFP}}^{\rm{\text{Mar}}}=\frac{\gamma}{k_{\rm{B}}T}\int\limits_{Q_{\text{u}}}^{Q_{\text{b}}}dx\text{e}^{U(x)/k_{\rm{B}}T}\int\limits_{-\infty}^{x}dy\text{e}^{-U(y)/k_{\rm{B}}T}.\end{split}$
(2)

Deviations between $\tau_{\rm{MFP}}^{\rm{Mar}}$ and
$\tau_{\rm{MFP}}^{\rm{MD}}$ can be due to either neglecting position dependent
friction, or the absence of memory effects. It was shown recently that a
position dependent Markovian friction cannot consistently describe both
folding and unfolding [17, 11], suggesting that the position-dependent
Markovian friction model is not suitable for protein folding (note that this
analysis was carried out for the neutral $\alpha$3D protein considered in this
article [11]). A comparison of $\tau_{\rm{MFP}}^{\rm{Mar}}$ and
$\tau_{\rm{MFP}}^{\rm{MD}}$ for the two pH systems shows that Eq. 2 accurately
predicts the transition time for the low pH system (Fig. 3C). This is expected
since the memory time scale approaches $1{\times}10^{-2}\tau_{\text{D}}$, so
$\alpha$3D exhibits effectively Markovian folding, satisfying a model with
position-independent, memoryless friction. Folding in the neutral pH system,
however, is accelerated by a factor of 6 compared to the Markovian prediction
(Fig. 3C).

Figure 4: Non-Markovian barrier crossing predictions over a range of memory
time scales. Predictions by an asymptotic cross-over formula
$\tau^{\rm{asy}}_{\alpha}/\tau^{\rm{asy}}_{\text{HF}}$ (solid lines) and the
Grote-Hynes theory $\tau^{\rm{GH}}_{\alpha}/\tau^{\rm{GH}}_{\text{HF}}$
(dashed lines) are parametrised by a four-component exponential fit for the
memory kernels in Fig. 2A. All memory times are multiplied by a common factor
$\alpha$ and each prediction is normalised by the respective high-friction
Markovian limit ($\tau^{\rm{asy}}_{\text{HF}}$ and
$\tau^{\rm{GH}}_{\text{HF}}$). The coloured symbols are for
$\tau_{\rm{MD}}^{\rm{MFP}}/\tau_{\rm{Mar}}^{\rm{MFP}}$, located at
$\alpha\tau_{\text{mem}}$ corresponding to $\alpha=1$.

Barrier-crossing acceleration is a hallmark of non-Markovian effects and is
predicted by various non-Markovian reaction rate theories. In Fig. 4, we
compare simulation results to various non-Markovian predictions, confirming
that the neutral pH system folds rapidly, despite its significantly higher
friction, due to memory speed-up effects. We parametrize both the Grote-Hynes
(GH) theory [32] and a recent asymptotic cross-over formula [10, 33, 34] by
fitting memory kernels to
$\Gamma(t)=\sum_{n=1}^{4}\gamma_{n}\text{e}^{-t/\tau_{n}}/\tau_{n}$
(Supplementary Materials Sec. IV) and uniformly scaling all memory times that
enter each prediction by a common factor $\alpha$.

$\tau^{\rm{asy}}_{\alpha}/\tau^{\rm{asy}}_{\text{HF}}$ is the asymptotic
cross-over formula [33, 34], rescaled by the high-friction, or overdamped,
Markovian limit $\tau^{\rm{asy}}_{\text{HF}}$ (Supplementary Material Sec. V).
This model was recently shown to be accurate for predicting the folding and
unfolding times of a large population of fast-folding proteins [11], even
capturing the non-Markovian barrier-crossing slow-down regime in some
instances. Similarly, $\tau^{\rm{GH}}_{\alpha}/\tau^{\rm{GH}}_{\text{HF}}$ is
the Grote-Hynes prediction, rescaled by the high-friction limit (Supplementary
Material Sec. V). In Fig. 4, we show both predictions, plotted as a function
of the 1st-moment memory time associated with the $\alpha$ scaling
($\alpha\tau_{\text{mem}}$). Both models predict reaction acceleration as
$\tau_{\text{mem}}/\tau_{\text{D}}>1{\times}10^{-2}$. By including the values
for $\tau_{\rm{MD}}^{\rm{MFP}}/\tau_{\rm{Mar}}^{\rm{MFP}}$ given in Fig. 3C,
we see that the deviations from Markovianity exhibited by the neutral pH
protein system are expected since this system is deep in the memory-induced
speed-up regime, as predicted by both non-Markovian models. Likewise, the low
pH system aligns well with the Markovian limit, with deviations being likely
due to the fact that both non-Markovian models incorporate harmonic
approximations to the free energy profiles, which are not well represented for
the two $\alpha$3D systems.

In summary, the extensive all-atom simulations of $\alpha$3D under two pH
conditions [24, 25] have provided us with an ideal data set to explore the
role of non-Markovian memory effects in protein folding dynamics. The
reduction in pH, which physiochemically results in a reduction in the inter-
residue salt-bridge interactions dispersed throughout the 73-residue chain,
manifests as a direct perturbation on the friction acting on the coarse-
grained reaction coordinate. By directly extracting the friction from the time
series of the $Q$ reaction coordinate, we have shown that this perturbation
modulates the amplitude and decay of the time-dependent friction kernel and
that, in this instance, the time scale is reduced by so much that the pH=2
system behaves as an effectively Markovian system. For the neutral pH system,
this is not the case and memory effects need to be explicitly accounted for to
predict the dramatic decrease in folding time compared to that predicted by
purely Markovian models.

We acknowledge support by the ERC Advanced Grant No. 835117 NoMaMemo and the
Deutsche Forschungsgemeinschaft (DFG) Grant No. SFB 1078. We thank the group
of David E. Shaw for providing the all-atom simulation data analysed in this
paper.

## References

  * Camacho and Thirumalai [1993] C. J. Camacho and D. Thirumalai, Proceedings of the National Academy of Sciences 90, 6369 (1993).
  * Bryngelson  _et al._ [1995] J. D. Bryngelson, J. N. Onuchic, N. D. Socci, and P. G. Wolynes, Proteins: Structure 21 (1995).
  * Dill and Chan [1997] K. A. Dill and H. S. Chan, Nature Structural Biology 4, 10 (1997).
  * Plotkin and Wolynes [2003] S. S. Plotkin and P. G. Wolynes, Proceedings of the National Academy of Sciences 100, 4417 (2003).
  * Sosnick and Barrick [2011] T. R. Sosnick and D. Barrick, Current Opinion in Structural Biology 21, 12 (2011).
  * Hinczewski  _et al._ [2013] M. Hinczewski, J. C. M. Gebhardt, M. Rief, and D. Thirumalai, Proceedings of the National Academy of Sciences 110, 4500 (2013).
  * Muñoz and Eaton [1999] V. Muñoz and W. A. Eaton, Proceedings of the National Academy of Sciences 96, 11311 (1999).
  * Best and Hummer [2006] R. B. Best and G. Hummer, Physical Review Letters 96, 228104 (2006).
  * Zheng and Best [2015] W. Zheng and R. B. Best, The Journal of Physical Chemistry B 119, 15247 (2015).
  * Kappler  _et al._ [2018] J. Kappler, J. O. Daldrop, F. N. Brünig, M. D. Boehle, and R. R. Netz, The Journal of Chemical Physics 148, 14903 (2018).
  * Dalton  _et al._ [2023a] B. A. Dalton, C. Ayaz, H. Kiefer, A. Klimek, L. Tepper, and R. R. Netz, Proceedings of the National Academy of Sciences 120, e2220068120 (2023a).
  * Plotkin and Wolynes [1998] S. S. Plotkin and P. G. Wolynes, Physical Review Letters 80, 5015 (1998).
  * Berezhkovskii and Makarov [2018] A. M. Berezhkovskii and D. E. Makarov, The Journal of Physical Chemistry Letters 9, 2190 (2018).
  * Berne and Harp [1970] B. J. Berne and G. D. Harp, On the Calculation of Time Correlation Functions (1970).
  * Daldrop  _et al._ [2018] J. O. Daldrop, J. Kappler, F. N. Brünig, and R. R. Netz, Proceedings of the National Academy of Sciences 115, 5169 (2018).
  * Kowalik  _et al._ [2019] B. Kowalik, J. O. Daldrop, J. Kappler, J. C. F. Schulz, A. Schlaich, and R. R. Netz, Physical Review E 100, 12126 (2019).
  * Ayaz  _et al._ [2021] C. Ayaz, L. Tepper, F. N. Brünig, J. Kappler, J. O. Daldrop, and R. R. Netz, Proceedings of the National Academy of Sciences 118, e2023856118 (2021).
  * Brünig  _et al._ [2022a] F. N. Brünig, O. Geburtig, A. von Canal, J. Kappler, and R. R. Netz, The Journal of Physical Chemistry B 126, 1579 (2022a).
  * Brünig  _et al._ [2022b] F. N. Brünig, J. O. Daldrop, and R. R. Netz, The Journal of Physical Chemistry B 126, 10295 (2022b).
  * Dalton  _et al._ [2023b] B. A. Dalton, H. Kiefer, and R. R. Netz, arXiv arXiv:2309.07521 (2023b).
  * Walsh  _et al._ [1999] S. T. R. Walsh, H. Cheng, J. W. Bryson, H. Roder, and W. F. DeGrado, Proceedings of the National Academy of Sciences 96, 5486 (1999).
  * Zhu  _et al._ [2003] Y. Zhu, D. O. V. Alonso, K. Maki, C.-Y. Huang, S. J. Lahr, V. Daggett, H. Roder, W. F. DeGrado, and F. Gai, Proceedings of the National Academy of Sciences 100, 15486 (2003).
  * Chung  _et al._ [2011] H. S. Chung, I. V. Gopich, K. McHale, T. Cellmer, J. M. Louis, and W. A. Eaton, The Journal of Physical Chemistry A 115, 3642 (2011).
  * Lindorff-Larsen  _et al._ [2011] K. Lindorff-Larsen, S. Piana, R. O. Dror, and D. E. Shaw, Science 334, 517 LP (2011).
  * Chung  _et al._ [2015] H. S. Chung, S. Piana-Agostinetti, D. E. Shaw, and W. A. Eaton, Science (New York, N.Y.) 349, 1504 (2015).
  * Best  _et al._ [2013] R. B. Best, G. Hummer, and W. A. Eaton, Proceedings of the National Academy of Sciences 110, 17874 LP (2013).
  * Hummer [2004] G. Hummer, The Journal of Chemical Physics 120, 516 (2004).
  * Shakhnovich  _et al._ [1991] E. Shakhnovich, G. Farztdinov, A. M. Gutin, and M. Karplus, Physical Review Letters 67, 1665 (1991).
  * Zwanzig [1961] R. Zwanzig, Physical Review 124, 983 (1961).
  * Mori [1965] H. Mori, Progress of Theoretical Physics 33, 423 (1965).
  * Zwanzig [2001] R. Zwanzig,  _Nonequilibrium Statistical Mechanics_ (Oxford University Press, 2001).
  * Grote and Hynes [1980] R. F. Grote and J. T. Hynes, The Journal of Chemical Physics 73, 2715 (1980).
  * Kappler  _et al._ [2019] J. Kappler, V. B. Hinrichsen, and R. R. Netz, The European Physical Journal E 42, 119 (2019).
  * Lavacchi  _et al._ [2020] L. Lavacchi, J. Kappler, and R. R. Netz, EPL 131 (2020).


# All the information in the integers is in the primes

Aidan Rocke
<EMAIL_ADDRESS>

###### Abstract

Building upon the work of Chebyshev, Shannon and Kontoyiannis, it may be
demonstrated that Chebyshev’s asymptotic result:

$\ln N\sim\sum_{p\leq N}\frac{1}{p}\cdot\ln p$

has a natural information-theoretic interpretation as the information-content
of a typical integer $Z\sim U([1,N])$, sampled under the maximum entropy
distribution, is asymptotically equivalent to the information content of all
primes $p\leq N$ weighted by their typical frequency.

##  1 Theorem: All the Statistical Information in the integers is contained in
the primes

The strategy in this section is to derive Chebyshev’s theorem(1852) via
Maximum Entropy inference, using an adaptation of Kontoyiannis’ original
arguments [4]. This theorem is an important precursor to the Prime Number
Theorem.

As the Unique Prime Factorisation theorem implies that the prime numbers
define a uniquely-decodable code for the integers, $Z\sim U([1,N])$
corresponds to a random prime factorisation in terms of the random variables
$X_{p}$:

$\forall Z\sim U([1,N]),Z=\prod_{p\leq N}p^{X_{p}}\implies H(Z)=\sum_{p\leq
N}H(X_{p})=\ln N$ (1)

Now, under Maximum Entropy assumptions the integers have a stochastic encoding
in terms of the prime numbers given by $P(Z\bmod p=0)\sim\frac{1}{p}$ since
for a given mean $\mathbb{E}[X_{p}]$, the geometric distribution maximises the
entropy of the integers:

$\mathbb{E}[X_{p}]=\sum_{k\geq 1}P(X_{p}\geq k)=\sum_{k\geq
1}\frac{1}{N}\lfloor\frac{N}{p^{k}}\rfloor\sim\frac{1}{p}$ (2)

so we have:

$H(Z)=-\ln\prod_{p\leq N}P(Z\bmod p=0)^{P(Z\bmod p=0)}\sim\sum_{p\leq
N}\frac{1}{p}\cdot\ln p$ (3)

By combining (1) and (3) we may deduce Chebyshev’s theorem:

$\sum_{p\leq N}\frac{1}{p}\cdot\ln p\sim\ln N$ (4)

which implies that almost all integers have a random prime factorisation and
that the average information gained from observing a prime number in the
interval $[1,N]$, weighted by their typical frequency, is on the order of
$\sim\ln N$.

It is worth noting that the information-theoretic formulation of Chebyshev’s
theorem that has been derived leads to a natural definition of the Statistical
Information of a prime number. This definition is consistent with the
hypothesis that the prime numbers have a uniform source distribution [1]. In
fact, the information-theoretic derivation of Chebyshev’s theorem and the
Prime Number theorem are consequences of this hypothesis.

##  2 Definition: The Statistical Information of a Prime Number

Given the information-theoretic derivation of the Prime Number Theorem,
presented in [1], the average Statistical Information in a prime number in the
interval $[1,N]$ is asymptotically on the order of $\sim\ln N$ or to be
precise,

$\frac{\log_{2}2^{N}}{\pi(N)}\sim\ln N$ (5)

as each integer in the interval $[1,N]$ is either prime or not prime.

##  3 Chebyshev’s theorem via Occam’s razor

Curiously, Chebyshev’s theorem may also be derived using Solomonoff’s theory
of Algorithmic Probability [10,11] which is conventionally known as Occam’s
razor.

As almost all integers are algorithmically random, $Z\sim U([1,N])$ might as
well be encoded using $\sim\log_{2}Z$ coin flips. It follows that the
algorithmic probability of observing a prime $p\in\mathbb{P}$ may be deduced
from Levin’s coding theorem:

$K_{U}(x)\sim-\log_{2}m(x)\implies m(p)\sim 2^{-K_{U}(p)}\sim
2^{-\log_{2}p}=\frac{1}{p}$ (6)

This has a natural frequentist interpretation, as we have a Bernoulli
distribution:

$\forall Z\sim U([1,N]),P(Z\bmod p=0)+P(Z\bmod p\neq
0)=\frac{1}{p}+\big{(}1-\frac{1}{p}\big{)}=1$ (7)

Now, using the fact that the Expected Kolmogorov Complexity of a random
variable equals its Shannon Entropy [5]:

$\mathbb{E}[K_{U}(Z)]=H(Z)+\mathcal{O}(1)$ (8)

we may deduce that:

$\forall Z\sim U([1,N]),\mathbb{E}[K_{U}(Z)]\sim-\log_{2}\prod_{p\leq
N}P(Z\bmod p=0)^{P(Z\bmod p=0)}\sim\sum_{p\leq
N}\frac{1}{p}\cdot\log_{2}p\sim\log_{2}N$ (9)

which yields Chebyshev’s theorem in base-2.

##  4 Discussion

In this article, we demonstrated that both Occam’s razor and Maximum Entropy
inference may be used to derive Chebyshev’s theorem for the distribution of
primes. The effective principle appears to be that the Expected Kolmogorov
Complexity of a random variable equals its Shannon Entropy.

There are good reasons to believe that this common principle may be developed
into a theory which may be applied to other non-trivial problems in
probabilistic number theory.

Acknowledgements: I would like to thank Anders Södergren for helpful
discussions regarding the structure of this publication, as well as Ioannis
Kontoyiannis for encouraging me to consider the relation between [1] and [4].

## References

[1] Aidan Rocke. AN INFORMATION-THEORETIC UPPER BOUND ON PRIME GAPS. Arxiv.
2021.

[2] P.L. Chebychev. Mémoire sur les nombres premiers. J. de Math. Pures Appl.,
17:366–390, 1852.

[3] P.L. Chebychev. Sur la totalité des nombres premiers inférieurs à une
limite donnée. J. de Math. Pures Appl., 17:341–365, 1852.

[4] I. Kontoyiannis. "Some information-theoretic computations related to the
distribution of prime numbers." In Festschrift in Honor of Jorma Rissanen, (P.
Grunwald, P. Myllymaki, I. Tabus, M. Weinberger, B. Yu, eds.), pp. 135-143,
Tampere University Press, May 2008.

[5] Peter Grünwald and Paul Vitányi. Shannon Information and Kolmogorov
Complexity. 2010.

[6] Lance Fortnow. Kolmogorov Complexity. 2000.

[7] Peter D. Grünwald. The Minimum Description Length Principle . MIT Press.
2007.

[8] A. N. Kolmogorov Three approaches to the quantitative definition of
information. Problems of Information and Transmission, 1(1):1–7, 1965

[9] G. J. Chaitin On the length of programs for computing finite binary
sequences: Statistical considerations. Journal of the ACM, 16(1):145–159,
1969.

[10] R. J. Solomonoff A formal theory of inductive inference: Parts 1 and 2.
Information and Control, 7:1–22 and 224–254, 1964.

[11] Marcus Hutter et al. (2007) Algorithmic probability. Scholarpedia,
2(8):2572.

[12] E.T. Jaynes. Information Theory and Statistical Mechanics. The Physical
Review. 1957.