\section{Full Details and Analysis of Theoretically Efficient {\textsc{LaB-PQ}}\xspace}\label{sec:app-tour-tree-cost}

We now give the proof of \cref{thm:main-tour-tree}.
We first analyze the functions in \cref{alg:tour-tree}.

\begin{lemma}
\label{lem:wtfunctions}
Consider a tournament tree\xspace{} on a universe of $n$ records, implemented with algorithms in \cref{alg:tour-tree}.
Given $b$ invocations of \mf{Mark}\xspace{} function in a batch, the total cost of these \mf{Mark}\xspace{} functions is $O\left(b\log\frac{n}{b}\right)$.
If there are $b$ invocations of \mf{Mark}\xspace{} function in the last batch, the total cost of the \mf{Sync}\xspace{} is $O\left(b\log\frac{n}{b}\right)$.
If the \mf{ExtractFrom}\xspace{} algorithm extracts $b$ (smallest) elements from tournament tree\xspace{} $T$, the total cost of \mf{ExtractFrom}\xspace{} is $O\left(b\log\frac{n}{b}\right)$.
\end{lemma}
\begin{proof}
  We first show that, in a tournament tree\xspace{}, for a subset $X$ of tree leaves, if we denote $S(X)$ as the set of all ancestors of nodes in $X$, then $|S(X)|=O\left(|X|\log\frac{n}{|X|}\right)$.
  This has been shown on more general self-balanced binary trees \cite{Blelloch2016justjoin,blelloch2018geometry}, and just a simplified case suffices as tournament tree\xspace{s} are complete binary trees.

  First, \mf{Mark}\xspace{} modifies the flag $\mb{renew}$ for each tree path node all modified leaves to the root. Let $X$ be all tree leaves corresponding to the invocations \mf{Mark}\xspace{} functions in this batch. Therefore $|X|=b$.
  By definition of $S(X)$, we know that only nodes in $X\cup S(X)$ are visited.
  Because of the {\texttt{test\_and\_set}}{} operation, a node $v\in S(X)$ recursively call \mf{Mark}\xspace{} on its parent $u$ if and only if $v$ successfully set the $\mb{renew}$ mark of $u$. This means that for any interior node $u\in S(X)$, this can happen only once. Therefore, every node in $S(X)$ is visited by the \mf{Mark}\xspace{} function at most once. This proves that the cost of all \mf{Mark}\xspace{} functions $O\left(b\log\frac{n}{b}\right)$.

  For \mf{Sync}\xspace{}, it restores all relevant interior tree nodes top-down. Let $X$ be all tree leaves corresponding to the invocations \mf{Mark}\xspace{} functions in the last batch (so that they need to be addressed in this \mf{Sync}\xspace{}). Therefore $|X|=b$.
  Note that the algorithm skips a subtree when the $\mb{renew}$ flag is \emph{false}{}. Therefore, the total number of visited nodes must be a subset of all nodes in $X\cup S(X)$ and their children. Given that each node has at most two children, this number can also be asymptotically bounded by $|X\cup S(X)|$, which are those marked \emph{true}{} in the $\mb{renew}$ flags. This gives the same bound as the total cost of \mf{Mark}\xspace{} functions in a batch.

  For \mf{Extract}\xspace{}, denote the leaves corresponding to all the output record\xspace{s} as set $X$ of size $b$.
  Note that we skip a subtree if its key (minimum key in its subtree) is larger than $\theta$. This means that we visit an interior node only if at least one of its descendants will be included in the output batch. Therefore, all visited nodes are also $X\cup S(X)$ and all their children.
  For all visited leaves, \mf{Extract}\xspace{} also calls \mf{Mark}\xspace{} to set $\mb{inQ}$ flag to 0 (\emph{false}). As proved above, the total cost of these \mf{Mark}\xspace{} functions is $O\left(b\log\frac{n}{b}\right)$.
  Therefore, the total cost of \mf{ExtractFrom}\xspace{} is also $O\left(b\log\frac{n}{b}\right)$ to extract a batch of size~$b$.
\end{proof}

With \cref{lem:wtfunctions} that shows the cost of each function in \cref{alg:tour-tree}, we can now formally prove \cref{thm:main-tour-tree}.

\medskip
\begin{proof}[Proof of \cref{thm:main-tour-tree}]
  Recall that the \emp{modification work} on a batch $B$ as all work paid to \mf{Update}\xspace{} all record\xspace{s} in $B$, as well as any later work (possibly done by a later \mf{Extract}\xspace{}) to actually apply the updates, and the \emp{extraction work} on a batch $B$ is all work paid to output the batch from the \mf{Extract}\xspace{} function, as well as any later work (possibly done by the next \mf{Extract}\xspace{}) to actually remove them from $Q$. In tournament tree\xspace{}, the modification work includes all \mf{Mark}\xspace{} operations called by the \mf{Update}\xspace{} functions in the batch, as well as later cost in the \mf{Sync}\xspace{} operation to restore keys of all the relevant interior nodes. From \cref{lem:wtfunctions}, the total cost is $O(b\log(\frac{n}{b}))$ for a batch of size $B$.
  The extraction work on a batch $B$ includes the work done by \mf{ExtractFrom}\xspace{} to get the output sequence, as well as to restore keys of all the relevant interior nodes in the next \mf{Sync}\xspace{}.

  For both \mf{Extract}\xspace{} and \mf{Update}\xspace{}, the span is no more than the height of the tree, which is $O(\log n)$.
\end{proof}

\section{Lower Bound for Batch-Update Work}\label{app:lower-bound}

We now show that in order to preserve the total ordering of all records, we need at least $O(\log (n+b))$ work per update, where $n$ is the number of records in the priority queue.
For a batch update of $b$ unordered records, there are ${b+n\choose b} \cdot b!$ total possible input cases.
In the comparison model, the lower bound for work of the batch update is:
 \begin{align*}
     \log_2\left({{b+n\choose b} \cdot b!}\right)=&\log_2\left({(b+n)!\over b!n!}\cdot b!\right)\\
     =&\log_2\left({\prod_{i=1}^{b}{(n+i)}}\right)\\=&\Omega(b\log (n+b))
 \end{align*}

Hence, to achieve better bounds, like in \cref{thm:main-tour-tree}, we cannot maintain the total ordering of all records.

\section{Finding the $\rho$-th Element in a List}\label{app:cpam}

We now discuss how to find the $\rho$-th element in a list, which is used in $\rho$-Stepping\xspace.
One solution is to pick $s=c(n/\rho+\log n)$ random samples where $c>1$ is a constant.
We can then sort the samples and pick the $(\rho s/n)$-th one.
By setting up the parameters correctly and using Chernoff Bound, we can show that with high probability, the selected output is within $(1-\epsilon)\rho$-th and $(1+\epsilon)\rho$-th element for some constant $0<\epsilon<1$.
We can check if the output is in this range by calling \mf{Extract}\xspace of the {\textsc{LaB-PQ}}\xspace{}, and if not, we can just retry until succeed.
We note that instead of finding the exact $\rho$-th element, using an approximate $\rho'$-th element within a constant factor of $\rho$ will not affect any bounds for $\rho$-Stepping\xspace.
In practice, we set $c=10$ and run the entire sampling algorithm sequentially, and this simple approach always gives a satisfactory output.
However, the work and span bounds are high probability bounds, instead of deterministic bounds.
Also, it only applies to the case when $\rho=\Omega(\sqrt{n})$.
Otherwise, the sorting cost will dominate.

We are aware of a data structure~\cite{cpam} (unpublished work) that can find the $\rho$-th element in a list efficiently.
More accurately, it takes $O(\rho)$ work and $O(\log n)$ span to find an approximate $\rho$-th element, and $O(\log (n/b))$ work per element in a batch insertion or a batch deletion of $b$ elements.
Basically, this is a blocked linked list data structure that looks like a binary search tree, but each leaf contains a block of $[\rho,3\rho]$ unsorted elements.
When the leaf grows too big or too small, we split or merge it with an amortized constant work per update and $O(\log n)$ span.
Hence, the cost per update is $O(\log (n/b))$ (to find the leaf node) plus a constant (amortized work for future leaf split or merge).
Since the first leaf contains the smallest $\rho$ to $3\rho$ record\xspace{s}, we can just pick the largest key among this leaf.
This gives us a key ranked within $\rho$ and $3\rho$, with the exception if we have fewer than $\rho$ records in total (in this case the largest is returned).

Using this data structure requires some changes since we have to explicitly generate the batch of updates (relaxations).
Theoretically, this is achievable as described below, but in practice the overhead can be significant so we do not use it.
This step is essentially Ligra's approach to generate the next frontier~\cite{ShunB2013}.
For each step, we allocate an array, and the size is the number of total neighbors of the vertices processed in this step.
Then we mark a flag for each successful relaxation, and after all relaxations are completed, we pack this array, and this is the batch of updates that applies to the next frontier (the priority queue in our case).
This step takes linear work and logarithmic span, so the cost is asymptotically bounded by the relaxation cost.

\section{Applying {\textsc{LaB-PQ}}\xspace{} to Shi-Spencer}
\label{app:shi-spencer}

Shi-Spencer algorithm~\cite{Shi1999} is a parallel SSSP algorithm with theoretical guarantee.
This algorithm is complicated, which harms its practicability.
We note that if we use our tournament-tree-based {\textsc{LaB-PQ}}\xspace{s} to replace the parallel search trees in Shi-Spencer algorithm, we can improve the work bound by up to a logarithmic factor. This algorithm needs preprocessing to shortcut each vertex to its $\rho$ nearest vertices\footnote{The notation used in the original paper is $k$-nearest vertices, we use $\rho$ here to be consistent with the results in our paper.}.

The main bottleneck in Shi-Spencer algorithm is to maintain a binary search tree for each vertex that keeps track of its $\rho$ nearest neighbors.
Once a vertex is settled, it will be removed from all the lists of their neighbors.
Such removals can happen for $O(m+n\rho)$ times in total.
The algorithm requires querying and extracting the $\rho$ nearest vertices of a vertex, which is used for $O(n)$ times in total.
In their original algorithm, they use parallel 2-3 trees to maintain the closest neighbor sets, which incur a logarithmic cost per update.

We note that the functionalities can be implemented by our tournament tree\xspace-based {\textsc{LaB-PQ}}\xspace{s}.
More importantly, similar to the stepping algorithms, since there are more updates than extractions, we can apply the updates lazily for better work efficiency.
In this case, the worst case is all the updates distribute evenly---each extract needs to first lazily update $O(m/n+\rho)$ previous updates and then apply the $O(n)$ extractions.
By the distribution lemma (\cref{lem:distribute}), the total work is:
$$O\left((m+n\rho)\log\frac{n\cdot n}{m+n\rho}\right)=O\left((m+n\rho)\log\frac{n^2}{m+n\rho}\right)$$
Hence, for dense graphs or if the algorithm picks a large $\rho$, the new work can be improved by up to a logarithm factor.
We note that a few other search trees and data structures are maintained in Shi-Spencer, but the total cost is  asymptotically bounded by the cost discussed in the above paragraph.

\section{Fully Dynamic {\textsc{LaB-PQ}}\xspace}\label{sec:fully-dynamic}

In \cref{sec:bdpq} we discuss and analyze the {\textsc{LaB-PQ}}\xspace's implementations assuming the universe of the records has a fixed size $n$, which is $|V|$ for SSSP.
This is good enough to derive the bounds in \cref{tab:bounds}.
However, we note that it is not hard to extend our algorithms to deal with the case without this assumption.

To do so, we need an explicit batch of updates, which can include insertions, deletions, and/or just key updates.
Our algorithm will in turn process the same type of operations.
For insertions (say $k$ in total), we can grow the tree size from $2n-1$ to $2(n+k)-1$, copy the leaf nodes from the range of $n$ to $n+k-1$ to the range of $2n$ to $2n+k-1$, insert the new keys to the range of $2n+k$ to $2(n+k)-1$ ($k$ in total), and update the corresponding tree paths.
The cost for $k$ insertions is $O(k+\log n)$.
For deletions (again, say $k$ in total), we can use the last tree leaves to fill in the holes for the deletions, and update the tree paths.
The cost for $k$ deletions is $O(k\log (n/k))$.
For updates, we can directly apply \cref{alg:tour-tree}.
All of them can be parallelized and have $O(\log n)$ span by a divide-and-conquer approach to update tree paths.

\section{Implementation Details for Resizable Hash Tables}\label{app:resize-ht}

As mentioned in \cref{sec:implementation}, we use a resizable hash table \variablename{next\_frontier} to maintain the frontier in our sparse version of stepping algorithm implementations.
Whenever a vertex $v$ is to be added to the next frontier, we pick a slot $i$ uniformly at random, and attempt to write $v$ to \variablename{next\_frontier[i]} using {\texttt{compare\_and\_swap}}{}, and use linear probing if there is collision or conflict due to concurrency. During the process, we keep an estimation of the size of \variablename{next\_frontier}, and resize accordingly. For frontier size $s$, the space used by the hash table is always $O(s)$.
We use two hash tables, alternating them to be the current and next frontiers.

It is worth noting that, even in resizing we do not explicitly allocate memory or move data. We will start with allocating an array of size $n$, and adjust tail pointer of the array to control the current hash table size. When the load factor of the current hash table reaches a pre-defined constant, we enlarge the hash table. We do this by moving the tail pointer to be twice the current size, and let all future scatters go to the last half of the hash table.
For each step, we always start from a pre-defined \texttt{MIN\_SIZE} of the hash table. Since we only insert to the hash table, we do not need to shrink the hash table. We show a pseudocode below to illustrate the process.

\begin{verbatim}
hash_table {
  vertex* table;
  int offset, tail, est_size;
  double SAMPLE_RATE;
  int MIN_SIZE;

  hash_table(int n) {
    table = new vertex[n];
    init(); }

  increment_size() {
    if (flip_coin(SAMPLE_RATE)) fetch_and_add(&est_size);}

  init() {offset = 0; tail = MIN_SIZE;}

  insert(vertex u) {
    slot = random_slot(u) + offset;
    while(!CAS(&table[slot], empty, u))
      slot = next_slot(slot);
    increment_size();
    if (size > est_size) {
      offset = tail;
      tail = tail*2; }
  }
};

hash_table* current_frontier, next_frontier;

parallel_for (v in current_frontier) {
  parallel_for (u is v's out-neighbor) {
    if (v relaxes u) {
      next_frontier.insert(u);
    } } }
swap(current_frontier, next_frontier);
next_frontier->init();
\end{verbatim}

\hide{
\section{$k$-$\rho$ property for different graphs}
\label{sec:krho}

To better understand the properties of different graphs, we show the $\rho$-$k_\rho$ curves in \cref{fig:k-rho}. We note that computing the exact value of $k_\rho$ is expensive (which makes Radius-Stepping\xspace{} impractical), and we estimate the value using 100 samples.
As mentioned, this indicates how much potential parallelism we can get on the original graph (i.e., without shortcuts). We show the value of $k_\rho$ when $\rho$ is $\log n$, $\sqrt{n}$, $n/\log n$, $n/10$, and $n$, respectively. Note that $k_n$ is the shortest path tree depth.

On all scale-free network{s}, we observe that $k_n$ is $O(\log n)$ (about $2\log n$). Interestingly, all scale-free network{s} are $(\log n,\sqrt{n})$-graph, which means that almost all vertices can reach their $\sqrt{n}$ nearest vertices by $\log n$ hops.
This is not surprising for scale-free network{s}, which have some ``hubs'' that are well connected to other vertices.
These vertices are easy to be reached by any source in a few steps.
Once reached, they can reach a lot of other vertices, which quickly accumulate $\sqrt{n}$ nearest neighbors for any source.
The road graphs have a different $\rho$-$k_\rho$ pattern.
On GE and USA, it takes more than $100$ hops to reach $\sqrt{n}$ nearest vertices, and $O(\sqrt{n})$ steps to reach $O(n)$ nearest vertices.
We believe the result is reasonable since road graphs are (almost) planar graphs.
}

\section{Experiment on Different Implementations for {\textsc{LaB-PQ}}\xspace{}}
\label{sec:exp:winningtree}
We also use a simple experiment to test the relative performance of array-based and tournament tree\xspace-based {\textsc{LaB-PQ}}\xspace{}.
We note that for \mf{Update}\xspace{}, our array-based implementation uses a hash-table-like data structure, which is only $O(1)$ work per insertion, while tournament tree\xspace{} needs up to logarithmic work.
For \mf{Extract}\xspace{}, the array-based implementation needs $O(n)$ time, while tournament tree\xspace{} needs $O(\rho\log(n/\rho))$ time to extract the top $\rho$ record\xspace{s}.
To estimate which one will give better performance in practice, we simulate the \mf{Extract}\xspace{} function and compare the performance for the two data structures.

\begin{figure}[th]
  \centering
  \includegraphics[width=\columnwidth]{figures-new/tree-array/tree-array.pdf}
  \caption{\textbf{Running time of extracting $\rho$ elements from tournament tree\xspace{}- and array- based implementation of {\textsc{LaB-PQ}}\xspace{}, respectively.}}\label{fig:tournament}
\end{figure}

To test the performance of \mf{Extract}\xspace{}, we vary the value of $\rho$. We first initialize the priority queue with $n=10^8$ record\xspace{s}, which is about the same order of magnitude as the number of vertices in the graphs we tested.  Our tournament tree\xspace{} is also implemented in flat arrays, avoiding expensive pointer-based structure.
We accumulate 10 runs of \mf{Extract}\xspace{}. The result is shown in \cref{fig:tournament}. The running time of array-based implementation is stable and is always around 0.35s. This is because the cost of array-based {\textsc{LaB-PQ}}\xspace{} is $O(n)$, which is not affected by the value of $\rho$. The cost of the tournament tree\xspace{} increases with the value of $\rho$.
When $\rho$ is larger than $2^{19}$, the array-based implementation achieves better performance.
We note that the setting is advantageous to tournament tree\xspace{}-based implementation since we do not consider the update cost ($O(1)$ for arrays and $O(\log n)$ for tournament tree\xspace{s}), and we always assume the frontier size to be large ($10^8$).
For the social networks we tested, the value of $\rho$ is usually larger than $2^{19}$. Therefore we always use array-based implementation of {\textsc{LaB-PQ}}\xspace{}. For the road graphs, the value of $\rho$ is slightly below $2^{19}$. As mentioned, this experiment does not consider update cost. Also, on road networks, the frontier size is also smaller (much smaller than $10^8$), which is favorable to array-based implementation. Based on these reasons, we choose to use the array-based data structures to implement our stepping algorithms.
From \cref{fig:visted-per-step}, for $\Delta$-Stepping\xspace{}, in the dense rounds (which dominate total running time), the number of vertices to be extracted is also reasonably large.
It is an interesting future work to see if tournament tree\xspace{}-based implementation can be more efficient on certain graphs that prefer a small $\rho$ in $\rho$-Stepping\xspace{}.

\section{Experiment on different source vertices}

In previous experiments, we show the average running time on multiple source vertices on each graph.
Now we show if different algorithms perform differently from the sources on the same graph.
To do so, we randomly pick 100 source vertices in Friendster and Twitter, and test the performance from each of them.
The results in shown in \cref{fig:source}.
Although the running time on these vertices and graphs differs, we can see that $\rho$-Stepping\xspace{} is consistently faster than the other two implementations, except for four vertices in Twitter.
We conclude that starting from fringe vertices and hub vertices does not make a significant difference in the relative performance between Bellman-Ford, $\Delta^*$-stepping, and $\rho$-Stepping\xspace{}.

\begin{figure}[ht]
  \centering
  \begin{tabular}{c}
    \includegraphics[width=\columnwidth]{figures-new/source/source-FT.pdf} \\
    \bf (a). FT  \\
    \includegraphics[width=\columnwidth]{figures-new/source/source-TW.pdf} \\
    \bf (b). TW \\
  \end{tabular}
  \caption{\textbf{Running time on difference source vertices.}}\label{fig:source}
\end{figure}

\section{Experiment on a different machine}

We also experiment on whether different machines affect the best $\rho$ choice. We launched a virtual machine in Amazon Web Services (AWS).
This single-socket machine is equipped with Intel Xeon Platinum 8175M CPUs with a total of 24 cores (48 hyperthreads). It has 192GB of main memory and 33MB L3 cache.
Other setting is the same as the previous experiments as mentioned in \cref{sec:exp:setup}. We present the result in \cref{fig:AWS-rho-time}.
Surprisingly, we found out that taking $\rho=2^{21}$ (or any value between $2^{20}$ to $2^{22}$) is still a good choice for all graphs for this different machine, although it has
much fewer cores.
Comparing \cref{fig:AWS-rho-time,fig:rho-time}, the machine with fewer cores shows a less stable performance especially when $\rho$ is small.
We note that there are also other parameters in the implementation that may affect the pattern of the relative performance of different $\rho$ values (e.g., the sparse-dense threshold).
In general, we do note that different hardware settings (e.g., cache size, number of sockets, number of threads) can affect the relative performance with different values of $\rho$, which is more significant for the smaller values of $\rho$.

\begin{figure}[ht]
  \centering
  \includegraphics[width=\columnwidth]{figures-new/AWS/AWS-rho-time.pdf}
  \caption{\textbf{Relative running time of $\rho$-Stepping\xspace{} with varied $\rho$ on the AWS machine.}}\label{fig:AWS-rho-time}
\end{figure}

\begin{figure*}[t]
  \centering
  \begin{tabular}{cccc}

    \multicolumn{4}{c}{\bf Friendster (FT):}\\
    \bf GAPBS &\bf Galois&\bf Julinne& \bf Ours (\textit{PQ-$\Delta^{*}$}\xspace)\\
    \includegraphics[width=0.5\columnwidth]{figures-new/delta/source/GAPBS-FT.pdf} &
    \includegraphics[width=0.5\columnwidth]{figures-new/delta/source/Galois-FT.pdf} &
    \includegraphics[width=0.5\columnwidth]{figures-new/delta/source/Julienne-FT.pdf} &
    \includegraphics[width=0.5\columnwidth]{figures-new/delta/source/ours-FT.pdf}
      \\
    \multicolumn{4}{c}{\bf WebGraph (WB):}\\
    \bf GAPBS &\bf Galois&\bf Julinne& \bf Ours (\textit{PQ-$\Delta^{*}$}\xspace)\\
    \includegraphics[width=0.5\columnwidth]{figures-new/delta/source/GAPBS-WB.pdf} &
    \includegraphics[width=0.5\columnwidth]{figures-new/delta/source/Galois-WB.pdf} &
    \includegraphics[width=0.5\columnwidth]{figures-new/delta/source/Julienne-WB.pdf} &
    \includegraphics[width=0.5\columnwidth]{figures-new/delta/source/ours-WB.pdf}
      \\
  \end{tabular}
  \caption{\textbf{Relative running time with varying $\Delta$ of multiple $\Delta$-Stepping\xspace{} implementations on different sources.} The five lines with different colors represent five sources. For each source, we normalize the running time to the corresponding fastest time. For each of the implementations, we use a window of $2^{11}$ around the best value of $\Delta$.  The best value of $\Delta$ is relatively stable for all implementations, except for a few instances (e.g., GAPBS and Julienne on WebGraph).
  Generally speaking, the best $\Delta$ for one source node makes another node at most 20\% slower.
  Also, it seems the result on directed graphs (WB) is more unstable than undirected graphs (FT).
  On both WB and FT, our \textit{PQ-$\Delta^{*}$}\xspace{} based on $\Delta^*$-Stepping\xspace{} shows very stable performance.

  }\label{fig:differentsource}
\end{figure*}

\begin{figure*}[t]
  \begin{tabular}{cccc}

    \multicolumn{4}{c}{\includegraphics[width=1.3\columnwidth]{figures-new/delta/caption.pdf}}\\
    \includegraphics[width=0.5\columnwidth]{figures-new/delta/delta_orkut.pdf} & \includegraphics[width=0.5\columnwidth]{figures-new/delta/delta_LJ.pdf} & \includegraphics[width=0.5\columnwidth]{figures-new/delta/delta_twitter.pdf} &
    \includegraphics[width=0.5\columnwidth]{figures-new/delta/delta_friendster.pdf} \\
\bf (a). OK& \bf (b). LJ & \bf (c). TW & \bf (d). FT \\
  \end{tabular}
  \begin{tabular}{ccc}
\includegraphics[width=0.5\columnwidth]{figures-new/delta/delta_web.pdf} & \includegraphics[width=0.5\columnwidth]{figures-new/delta/delta_GE.pdf} &
\includegraphics[width=0.5\columnwidth]{figures-new/delta/delta_US.pdf}\\
\bf (e). WB&\bf (f). GE&\bf (g). USA\\
  \end{tabular}
  \caption{\textbf{$\Delta$-stepping relative running time with varying $\Delta$.} We use 96 cores (192 hyperthreads).}\label{fig:alldelta}

\end{figure*}

\begin{figure*}
\begin{tabular}{cccc}
  \includegraphics[width=0.5\columnwidth]{figures-new/per-round/OK.pdf}&
  \includegraphics[width=0.5\columnwidth]{figures-new/per-round/LJ.pdf}&
  \includegraphics[width=0.5\columnwidth]{figures-new/per-round/TW.pdf}&
  \includegraphics[width=0.5\columnwidth]{figures-new/per-round/FT.pdf}\\
  \includegraphics[width=0.5\columnwidth]{figures-new/per-round/OK-e.pdf}&
  \includegraphics[width=0.5\columnwidth]{figures-new/per-round/LJ-e.pdf}&
  \includegraphics[width=0.5\columnwidth]{figures-new/per-round/TW-e.pdf}&
  \includegraphics[width=0.5\columnwidth]{figures-new/per-round/FT-e.pdf}\\
\bf (a). OK& \bf (b). LJ & \bf (c). TW& \bf (d). FT \\
\end{tabular}
\begin{tabular}{ccc}
  \includegraphics[width=0.5\columnwidth]{figures-new/per-round/WB.pdf}&
  \includegraphics[width=0.5\columnwidth]{figures-new/per-round/USA.pdf}&
  \includegraphics[width=0.5\columnwidth]{figures-new/per-round/GE.pdf}\\
  \includegraphics[width=0.5\columnwidth]{figures-new/per-round/WB-e.pdf}&
  \includegraphics[width=0.5\columnwidth]{figures-new/per-round/USA-e.pdf}&
  \includegraphics[width=0.5\columnwidth]{figures-new/per-round/GE-e.pdf}\\
\bf (e). WB& \bf (f). USA & \bf (g). GE\\
\end{tabular}
\caption{\textbf{Number of visited vertices in each step in \textit{PQ-$\rho$}\xspace, \textit{PQ-$\Delta^{*}$}\xspace and \textit{PQ-BF}\xspace.} \mdseries Here we only run on one source vertex, since it has unclear meaning to compute the average of multiple runs on each step.  Hence, the runtimes can be different from Table \ref{tab:alltime} (average on 100 runs from 10 source vertices), and some curves are bumpy. We use 96 cores (192 hyperthreads).\label{fig:visted-per-step-full}}
\end{figure*}

\section{{\textsc{LaB-PQ}}\xspace Implementation} \label{sec:bdpq}

We now discuss how to efficiently support {\textsc{LaB-PQ}}\xspace{} in \cref{alg:framework}.
We present two data structures for {\textsc{LaB-PQ}}\xspace with the goal of theoretical and practical efficiency, respectively.
The obliviousness for data structures from the algorithm's perspective is an advantage of the {\textsc{LaB-PQ}}\xspace ADT.

In our analysis,
we define a \emp{batch of modifications} as all \mf{Update}\xspace{} operations between two invocations of \mf{Extract}\xspace{} functions.
The \emp{modification work} on a batch $B$ is all work paid to \mf{Update}\xspace{} all record\xspace{s} in $B$, as well as any later work (done by a later \mf{Extract}\xspace{}) to actually apply the updates.
We define a \emp{batch of extraction} as all record\xspace{s} returned by an \mf{Extract}\xspace{} function.
The \emp{extraction work} on a batch $B$ is all work paid to output the batch from the \mf{Extract}\xspace{} function, as well as any later work (done by the next \mf{Extract}\xspace{}) to actually remove them from $Q$.

\subsection{Related Work}

Early PRAM and BSP algorithms had explored parallel priority queues in a variety of approaches~\cite{brodal1998parallel,chen1994fast,das1996optimal,pinotti1991parallel,ranade1994parallelism,crupi1996parallel,baumker1996realistic}, and heavily rely on synchronization-based techniques such as pipelining.
These algorithms do not have better bounds than recent batch-dynamic search trees~\cite{Blelloch2016justjoin,sun2018pam,blelloch2020optimal,sun2019implementing,sun2019parallel} when mapping to the fork-join model.
Other previous papers considered the concurrent, external-memory, and other settings~\cite{alistarh2015spraylist,sundell2005fast,linden2013skiplist,shavit2000skiplist,liu2012lock,henzinger2013quantitative,zhou2019practical,calciu2014adaptive,bingmann2015bulk,sanders2019sequential,sanders1998randomized,sanders2000fast}.
These data structures also do not have better bounds than batch-dynamic search trees since they do not focus on optimizing work or span.
However, existing batch-dynamic search trees or other data structures (e.g., skiplists) maintaining the total ordering of the record\xspace{s}, incur an $\Omega(\log(n))$ work lower bound per record\xspace{} update (more details are in the full paper).
Our key observation is that maintaining total ordering, which incurs overhead both theoretically and practically, is not necessary for a parallel priority queue.

To the best of our knowledge, the only parallel data structure that has similar bounds to our new data structure is the batch-dynamic binary heap~\cite{wang2020parallel}.
However, it has a few disadvantages: it does not support efficient batch-extract, is very complicated (no implementation available), and the span is suboptimal ($O(\log^2 n)$ in the binary fork-join model).
Our new tournament-tree based {\textsc{LaB-PQ}}\xspace{} supports full features in the {\textsc{LaB-PQ}}\xspace, has $O(\log n)$ span, and is arguably much simpler.

\subsection{Tournament-Tree-Based Implementation}\label{sec:tour-tree}

We start with introducing the tournament tree (aka.\ winner tree).
It is a complete binary tree with $n$ external nodes (leaves) and $n-1$ interior nodes.
A tournament tree\xspace stores the record\xspace{s} in the leaves.
In our use case, we only need to store the record\xspace{} ${\mathit{id}}\xspace$ in the leaves using the {\textsc{LaB-PQ}}\xspace interface.
Each interior node stores $k\in K$ ($K$ is in key type for the records) that takes the smaller key (defined by $<_K$) from its children.
\cref{fig:tour-tree}(a) illustrates a tournament tree\xspace when keys are integers and $<_K$ is $<_{\mathbb{Z}}$.

\input{tour-tree-illu}

We now discuss how to use a tournament tree\xspace to implement {\textsc{LaB-PQ}}\xspace{}.
We will use $t.\mb{left}$, $t.\mb{right}$ and $t.\mb{parent}$ to denote the left child, right child and parent of a node~$t$.
For simplicity, we assume the universe of the records has a fixed size $n$ (for SSSP $n=|V|$), and the tournament tree\xspace has $n$ leaf nodes each with a boolean flag $\mb{inQ}$ indicating if this record is in (has been inserted to) the {\textsc{LaB-PQ}}\xspace $Q$.
We note that this is sufficient for the SSSP algorithms.
The dynamic version (where the size of the tournament tree\xspace{} changes with the size of {\textsc{LaB-PQ}}\xspace{}) will be described in the full paper.

\begin{algorithm
\caption{The Tournament-tree based {\textsc{LaB-PQ}}\xspace.\label{alg:tour-tree}}
\fontsize{9pt}{10pt}\selectfont
\SetKw{KwMaintain}{Maintains}
\KwMaintain{A tournament tree\xspace $T$ with $n$ leaf nodes each corresponding to a record.}
    \vspace{0.5em}
\DontPrintSemicolon
\SetKwBlock{InPar}{In Parallel:}{}
\SetKwProg{myfunc}{Function}{}{}

\myfunc{\upshape $\mf{Mark}\xspace($record id ${\mathit{id}}\xspace$, boolean flag $\mb{newflag})$} {
  Let $t$ be the tree leaf corresponding to \mb{id}\\
  $t.\mathit{inQ}\xspace\gets \mb{newflag}$\\
  \While {\upshape $t\ne T.root$ and $\TAS(t.\mb{parent}.\mb{renew})$} {
    $t\gets t.\mb{parent}$
  }
}
{\vspace{.05in}}{}
\myfunc{\upshape $\mf{Sync}\xspace($node $t) \to k\in K$\label{line:update-start}} {
  \If {$t$ is leaf} {
    \lIf{$t.\mb{inQ}$} {\Return{$\delta[t.id]$}}
    \lElse {\Return{$+\infty$}
  }
  }
  \lIf {$t.\mb{renew}=0$} {\Return{$t.k$}
  $t.\mb{renew} \gets 0$\\
  \InPar {
    $\mb{leftKey}\gets \mf{Sync}\xspace(t.\mb{left})$\\
    $\mb{rightKey}\gets \mf{Sync}\xspace(t.\mb{right})$
  }
  \Return {$t.k \gets\min(\mb{leftKey}, \mb{rightKey}$)}\label{line:update-end}
}
{\vspace{.05in}}{}
\myfunc{\upshape $\mf{ExtractFrom}\xspace(\text{threshold }\theta\text{, node }t) \to \mb{seq}$\label{line:traverse-start}} {
  \If {$t$ is a leaf} {

    \If {$(\delta[t.{\mathit{id}}\xspace] \le \theta)$}{
    $\mf{Mark}\xspace(t.{\mathit{id}}\xspace, 0)$\label{line:delete} \tcp*{Marked as not in Q}
    \Return{$\{t.\mb{id}\}$}}
    \lElse{\Return{$\{\}$}}

  }
  \lIf {$\theta<t.k$} {\Return{$\{\}$}\tcp*[f]{empty seq}}
  \InPar {
    $\mb{leftseq}\gets\mf{ExtractFrom}\xspace(\theta,t.\mb{left})$\\
    $\mb{rightseq}\gets\mf{ExtractFrom}\xspace(\theta,t.\mb{right})$
  }
  \Return $\mb{leftseq}+\mb{rightseq}$\label{line:traverse-end}
}
{\vspace{.05in}}{}
\myfunc{\upshape $\mf{Extract}\xspace(\text{threshold }\theta$)} {
  $\mf{Sync}\xspace(T.\mb{root})$\\
  \Return $\mf{ExtractFrom}\xspace(\theta,T.\mb{root})$
}
{\vspace{.05in}}{}
\myfunc{\upshape $\mf{Update}\xspace({\mathit{id}}\xspace$)} {
  $\mf{Mark}\xspace({\mathit{id}}\xspace{},1)$\\
}
\end{algorithm}

A tournament tree\xspace $T$ on $n$ records contains $2n-1$ nodes in total. The first $n-1$ nodes are interior nodes.
For an interior node $t$, we use $t.k$ to denote the key stored in node $t$.
To support the {\textsc{LaB-PQ}}\xspace{} interface, each interior node contains a bit flag $\mb{renew}$ indicating if any key in its subtree has been modified after the last update of $t.k$.
This flag is initially set to 0 (\emph{false}).

Constructing such a tree with given initial values simply takes linear work and $O(\log n)$ span using divide-and-conquer: construct both subtrees recursively in parallel, and update the root's key based on the two children's keys.

We next present the implementation of {\textsc{LaB-PQ}}\xspace's interface using tournament tree\xspace{}. Due to page limit, we only show the pseudocode (\cref{alg:tour-tree}) and a high-level overview here.
The analysis are given in the full paper.
We first introduce a helper function \mf{Mark}\xspace{}$({\mathit{id}}\xspace{}, \mathtext{newflag})$.

\mf{\textbf{Mark}}\xspace{}$({\mathit{id}}\xspace{}, \mathtext{newflag})$ first sets the record ${\mathit{id}}\xspace$'s $\mathtext{inQ}$ flag to be $\mathtext{newflag}$---0 means the record should be deleted, and 1 means inserted.
Whichever value $\mathtext{newflag}$ is, this means the record\xspace{} of ${\mathit{id}}\xspace$ has been updated.
Then, the algorithm marks the $\mb{renew}$ flags of the nodes on the tree path from the updated leaf to the root.
This process is executed using \TAS{}.
If the \TAS{} fails, we know that another \mf{Mark}\xspace{} has marked the rest of the path, so the current \mf{Mark}\xspace{} terminates immediately.
An example is shown in \cref{fig:tour-tree}(b).
Updating the nodes' keys is postponed to the next \mf{Extract}\xspace{} function.

\para{\upshape \mf{\textbf{Update}}{}.}
The \mf{Update}\xspace{} algorithm simply calls \mf{Mark}\xspace{}$({\mathit{id}}\xspace{}, 1)$.

\para{\upshape \textbf{\textsc{Extract}}\xspace{}.} \mf{Extract}\xspace{} first uses a function \mf{Sync}\xspace{} to update the keys for all nodes with $\mb{renew}$ flag as 1.
It then calls \mf{ExtractFrom}\xspace{} to output all record\xspace{s} with keys no more than $\theta$.
Those output keys are also marked as deleted from $T$ (\cref{line:delete}).

The $\mf{Sync}\xspace{}(t)$ function recursively restores the keys in the interior nodes using a divide-and-conquer approach (\cref{line:update-start}--\ref{line:update-end}), and returns the key at the current node $t$.
The return value of a leaf node is either the record's key or infinity, depending on the $\mathit{inQ}\xspace$ flag.
For an interior node, \mf{Sync}\xspace{} will update the key to be the smaller one of its two children and return this key.
After all interior tree nodes have been updated, we use \mf{ExtractFrom}\xspace{} to acquire all records with keys no more than $\theta$.
This step can be parallelized similarly using divide-and-conquer (\cref{line:traverse-start}--\ref{line:traverse-end}): we can traverse the left and right subtrees respectively and concatenate the two results. A subtree is skipped when the key at the subtree root (minimum key in the subtree) is larger than $\theta$.
Note that if we want the output sequence in a consecutive array, we can traverse for two rounds---the first round computes the number of extracted records, and the second round writes them to the corresponding slots.

We can implement \mf{Reduce}\xspace similarly.
We keep a collective status $t.a\in A$ ($A$ is the augmented value type) for each interior node $t$, and it is updated in the \mf{Update} function in \cref{line:update-start}--\ref{line:update-end} similar to the update for $k$.
We do not need to update $t.a$ for node $t$ if the subtree rooted at $t$ remains unchanged, which is captured by $t.\mb{renew}$.

\begin{theorem}\label{thm:main-tour-tree}
  Consider a tournament tree\xspace{} on a universe of $n$ records, implemented with algorithms in \cref{alg:tour-tree}.
  The modification work on a size-$b$ batch is $O\left(b\log(n/b)\right)$.
  The extraction work on a size-$b$ batch is $O\left(b\log(n/b)\right)$.
  The span of \mf{Extract}\xspace{} and \mf{Update}\xspace{} is $O(\log n)$.
\end{theorem}

\ifx0\undefined
The formal analysis of \cref{thm:main-tour-tree} is given in
the full paper.
We will see how to use \cref{thm:main-tour-tree} in \cref{sec:sssp}.
\else
To prove this theorem, we will use the following lemma.
\begin{lemma}
\label{lem:wtfunctions}
Given $b$ invocations of \mf{Mark}\xspace{} function in a batch, the total cost of these \mf{Mark}\xspace{} functions is $O\left(b\log\frac{n}{b}\right)$.
If there are $b$ invocations of \mf{Mark}\xspace{} function in the last batch, the total cost of the \mf{Sync}\xspace{} is $O\left(b\log\frac{n}{b}\right)$.
If the \mf{ExtractFrom}\xspace{} algorithm extracts $b$ (smallest) elements from tournament tree\xspace{} $T$, the total cost of \mf{ExtractFrom}\xspace{} is $O\left(b\log\frac{n}{b}\right)$.
\end{lemma}
\begin{proof}
  We first show that, in a tournament tree\xspace{}, for a subset $X$ of tree leaves, if we denote $S(X)$ as the set of all ancestors of nodes in $X$, then $|S(X)|=O\left(|X|\log\frac{n}{|X|}\right)$.
  This has been shown on more general self-balanced binary trees \cite{Blelloch2016justjoin,blelloch2018geometry}, and just a simplified case suffices as tournament tree\xspace{s} are complete binary trees.

  First, \mf{Mark}\xspace{} modifies the flag $\mb{renew}$ for each tree path node all modified leaves to the root. Let $X$ be all tree leaves corresponding to the invocations \mf{Mark}\xspace{} functions in this batch. Therefore $|X|=b$.
  By definition of $S(X)$, we know that only nodes in $X\cup S(X)$ are visited.
  Because of the {\texttt{test\_and\_set}}{} operation, a node $v\in S(X)$ recursively call \mf{Mark}\xspace{} on its parent $u$ if and only if $v$ successfully set the $\mb{renew}$ mark of $u$. This means that for any interior node $u\in S(X)$, this can happen only once. Therefore, every node in $S(X)$ is visited by the \mf{Mark}\xspace{} function at most once. This proves that the cost of all \mf{Mark}\xspace{} functions $O\left(b\log\frac{n}{b}\right)$.

  For \mf{Sync}\xspace{}, it restores all relevant interior tree nodes top-down. Let $X$ be all tree leaves corresponding to the invocations \mf{Mark}\xspace{} functions in the last batch (so that they need to be addressed in this \mf{Sync}\xspace{}). Therefore $|X|=b$.
  Note that the algorithm skips a subtree when the $\mb{renew}$ flag is \emph{false}{}. Therefore, the total number of visited nodes must be a subset of all nodes in $X\cup S(X)$ and their children. Given that each node has at most two children, this number can also be asymptotically bounded by $|X\cup S(X)|$, which are those marked \emph{true}{} in the $\mb{renew}$ flags. This gives the same bound as the total cost of \mf{Mark}\xspace{} functions in a batch.

  For \mf{Extract}\xspace{}, denote the leaves corresponding to all the output record\xspace{s} as set $X$ of size $b$.
  Note that we skip a subtree if its key (minimum key in its subtree) is larger than $\theta$. This means that we visit an interior node only if at least one of its descendants will be included in the output batch. Therefore, all visited nodes are also $X\cup S(X)$ and all their children.
  For all visited leaves, \mf{Extract}\xspace{} also calls \mf{Mark}\xspace{} to set $\mb{inQ}$ flag to 0 (\emph{false}). As proved above, the total cost of these \mf{Mark}\xspace{} functions is $O\left(b\log\frac{n}{b}\right)$.
  Therefore, the total cost of \mf{ExtractFrom}\xspace{} is also $O\left(b\log\frac{n}{b}\right)$ to extract a batch of size~$b$.
\end{proof}

With \cref{lem:wtfunctions} that shows the cost of each function in \cref{alg:tour-tree}, we can now formally prove \cref{thm:main-tour-tree}.

\medskip
\begin{proof}[Proof of \cref{thm:main-tour-tree}]
  Recall that the \emp{modification work} on a batch $B$ as all work paid to \mf{Update}\xspace{} all record\xspace{s} in $B$, as well as any later work (possibly done by a later \mf{Extract}\xspace{}) to actually apply the updates, and the \emp{extraction work} on a batch $B$ is all work paid to output the batch from the \mf{Extract}\xspace{} function, as well as any later work (possibly done by the next \mf{Extract}\xspace{}) to actually remove them from $Q$. In tournament tree\xspace{}, the modification work includes all \mf{Mark}\xspace{} operations called by the \mf{Update}\xspace{} functions in the batch, as well as later cost in the \mf{Sync}\xspace{} operation to restore keys of all the relevant interior nodes. From \cref{lem:wtfunctions}, the total cost is $O(b\log(\frac{n}{b}))$ for a batch of size $B$.
  The extraction work on a batch $B$ includes the work done by \mf{ExtractFrom}\xspace{} to get the output sequence, as well as to restore keys of all the relevant interior nodes in the next \mf{Sync}\xspace{}.

  For both \mf{Extract}\xspace{} and \mf{Update}\xspace{}, the span is no more than the height of the tree, which is $O(\log n)$.
\end{proof}
\fi

\subsection{Array-Based Implementation}

\cref{alg:tour-tree} uses a tree-based structure to provide tight work bounds for applying a batch of modifications or extractions.
This is asymptotically better than batch-dynamic search trees~\cite{Blelloch2016justjoin,sun2018pam,blelloch2020optimal}.
However, in practice, maintaining a tree-based data structure can be expensive because of larger memory footprint and random access.
Even though we can implement a tournament tree\xspace{} in a flat array (no pointers), it still requires extra storage for interior nodes and incurs frequent random accesses (following tree path).
When the batch size $b$ approaches $n$ and $O(\log(n/b))$ becomes small, the theoretical advantage of tournament tree\xspace{s} becomes insignificant, and is asymptotically the same as just loop over all record\xspace{s}.

This is observed by the practitioners. Most (if not all) practical SSSP implementations just keep an array for all records without maintaining sophisticated structures.
This is because parallel SSSP algorithms usually use a very large value of $b$ to get sufficient parallelism.
To implement \mf{Update}\xspace{} on an array, we can just set a flag to indicate a record\xspace{} is added to $Q$.
For \mf{Extract}\xspace, we loop over the entire array and pack all records with keys within~$\theta$ in parallel, which takes linear work and $O(\log n)$ span.
While efficiently implementing the array requires many subtle details (shown in \cref{sec:impl-bdpq}), asymptotically, the following bound is easy to see.

\begin{theorem}\label{thm:array-based}
  The array-based {\textsc{LaB-PQ}}\xspace requires $O(b)$ modification work on a size-$b$ batch.
  The extraction work on a size-$b$ batch is $O(n)$.
  The span of \mf{Extract}\xspace{} and \mf{Update}\xspace{} is $O(\log n)$.
\end{theorem}

\section{Parallel Soft Priority Queue}
\label{sec:sssp}
Another typical example of abstract data type (ADT) is the priority queue (\pq{}).
A priority queue maintains a set of elements with keys as priorities, and supports three primitives: \mf{Extract-Min} that pops the element with the smallest key, \mf{Insert} that adds a new element, and \mf{Dec-Key} that decreases the key of an element.
Other operations such as \mf{Delete} or \mf{Update} can be implemented via these primitives.
Priority queues have many applications, such as in Prim's algorithm~\cite{prim1957shortest} in computing minimum spanning tree, constructing Huffman tree~\cite{huffman1952method}, and best-first searches~\cite{pearl1984intelligent}.
The most famous application of priority queues is probably in Dijkstra's algorithm that efficiently computes the single-source shortest path (SSSP\xspace{}) on a graph.
Sequentially, Dijkstra's algorithm is simple and practical.
Similar to the previous section, we will use SSSP and Dijkstra as a motivating example for parallel priority queues.

As oppose to many graph algorithms, SSSP is notoriously hard in the parallel setting.
To date, there is no work-efficient (or even close to) solution with non-trivial ($o(n)$) span bound.
Directly parallelizing Dijkstra's leads to $\Omega(n)$ span due to the use of priority queue: \mf{Extract-Min} pops one vertex at a time for $n$ times in total.
There are two major attempts to parallelize the priority queue.
One is to allow concurrent updates to the priority queue \cite{}.
Hence, processors can work on different vertices together in Dijkstra, but the \mf{Extract-Min}s from the priority queue are still sequentialized.
The other attempt is to relax the priority queue to allow batch updates and batch extracts on $k$ elements, so multiple processors can work on different vertices simultaneously.
To support these functions, the underlying data structure is a search tree, which can be either the parallel search trees by the PI and coauthors~\cite{Blelloch2016just,sun2018pam,sun2019parallel} in the shared-memory setting or B-trees \cite{} in the distributed setting.
Unfortunately, maintaining the search tree is very costly because the total ordering of all elements are kept, which is not required in the sequential heap implementations.

One could also use the parallel BST \cite{} proposed by one of the PIs as a batched priority queue (as is used in \cite{}). However, maintaining an extra tree structure is less I/O-friendly in practice.
Interestingly, the work-inefficient Bellman-Ford algorithm is easier to be parallelized: in each round, all vertices attempt to relax all their neighbors, which can be done fully in parallel. Previous work has shown that this algorithm is easy to implement and have reasonable performance in many cases \cite{}.
Although Bellman-Ford is far from work-efficient and is almost trivial to parallelize, to the date, there is no parallel SSSP\xspace{} algorithm that is more efficient than Bellman-Ford both in theory (better worst-case bounds) and practice (better performance).

Due to the problems in aforementioned parallel priority queues, both theoretical and practical parallel SSSP algorithms \emph{do not} use the priority queue abstraction.
In practice, Bellman-Ford can easily be parallelized.
In each {step}{}, only a subset of vertices are used to relax their neighbors, noted as the \emph{frontier}.
The algorithm maintains a queue buffering all vertices that recently have their distance updated, but have not yet relaxed their neighbors.
The frontier is extracted from the queue by a predefined condition (noted as $\mathcal{I}_{{\mathtext{ex}}}$), which usually selects vertices close to the source. After that, relaxed vertices will update their priority and be added to the queue.
In other words, these algorithms already use a conceptual \emph{Batch-PQ\xspace{}}, although not explicitly.
For example, the $\Delta$-Stepping\xspace~\cite{meyer2003delta} algorithm uses all vertices with tentative distance $\le i\times \Delta$ in the frontier for {step}{} $i$, where $\Delta$ is a graph-specific parameter. Due to its simplicity, $\Delta$-Stepping\xspace{} is the most widely-used parallel SSSP\xspace{} algorithm in practice, and is adopted in many graph processing systems \cite{dhulipala2017julienne,hang2020optimizing,nguyen2013lightweight,BeamerAP15}. $\Delta$-Stepping\xspace{} does not have worst-case time bound because the parameter $\Delta$ highly depends on the input graph structure and edge weight, and its best choice is not known before running the algorithm itself.
On the theory side, the state-of-the-art parallel SSSP\xspace{} theoretical result is the Radius-Stepping\xspace{} algorithm, which was proposed by the PIs.
Instead of pushing out the frontier by a fixed distance as in $\Delta$-Stepping\xspace{}, Radius-Stepping\xspace{} changes the parameter $\Delta$ dynamically, noted as the current \emph{radius} $r_i$. The radius is decided based on the current frontier. Using this idea, Radius-Stepping\xspace{} achieves worst-case work and span bounds and good parallelism (work-span ratio). There also exists other theory work using similar idea such as Shi and Spencers \cite{}. However, these algorithms are not used in practice because they require a preprocessing phase to add shortcuts. In practice, this could bring up serious space and time overhead.

We summarize the framework of frontier-based SSSP\xspace{} algorithm in Algorithm \ref{algo:sssp}.
All these high-level idea of these algorithms for SSSP\xspace{} shows similarity of employing a ``batched \pq{}'' conceptually (the \fname{SoftPQ}) structure in Algorithm \ref{algo:sssp}).
Namely, the condition function $\mathcal{I}_{{\mathtext{ex}}}$ defines the priority, and the frontier is selected based on this priority.
Instead of one element, $\mathcal{F}$ is extracted as batch.
$\mathcal{I}_{{\mathtext{ex}}}$ can be defined not only based on the value of a single key, but also more complex criteria. As such, we call this abstract data type a \emph{Soft Batched Priority Queue} (\emph{{Soft Batch-PQ}{}}).
Surprisingly, however, none of the previous work formally studies this abstraction and semantics of the parallel \pq{} for the SSSP\xspace{} problem.
This is probably due to the use of this batched \pq{} is not as straightforward: instead of extracting a batch and process it in every round (noted as {step}{s}) and repeat, there are also some ``{substep}{s}'' (Line \ref{line:substep}). When processing the current frontier, a relaxed vertex may be ``pulled'' into the frontier and will be processed in the same {step}{}, if some condition (noted as $\mathcal{I}_{{\mathtext{add}}}$) is satisfied. Usually this condition checks if the tentative distance is smaller enough. This allows the frontier to grow dynamically, and aims at reducing both work and span of the algorithm. First, the vertices with small tentative distance should be used to relax other vertices as early as possible to avoid redundant relaxation. Secondly, processing more vertices in the same round reduces the number of {step}{s}, and allows more work to do in parallel in the same {step}{}, which improves parallelism.

We argue that, this mechanism of allowing for the growth of the frontier, is the main feature that makes a priority queue ``batched'' instead of ``streaming'' (dealing with one element at a time), and ``parallel'' instead of ``serial'', in the SSSP\xspace{} problem.
In fact, most of the previous SSSP\xspace{} algorithms falls into this framework in Algorithm \ref{algo:sssp}. In Table \ref{tab:allsssp}, we list a few examples including the parallel Dijkstra's, parallel Bellman-Ford, $\Delta$-Stepping\xspace{}, Radius-Stepping\xspace{} and Shi and Spencer's algorithm. Their use of priority queue boils down to different choices of the two conditions: $\mathcal{I}_{{\mathtext{ex}}}$, which is the condition to extract the frontier, and $\mathcal{I}_{{\mathtext{add}}}$, which decides if a vertex should be added to the current frontier. We note that the pseudocode only presents the conceptual idea of these algorithms, which may not reflect the real implementation details.

In this project, by better understanding the semantics and designing an efficient implementation of the {Soft Batch-PQ}{}, our goal is to achieve an efficient parallel SSSP\xspace{} algorithm. We will first need to understand the challenge in the parallel SSSP\xspace{} algorithm, which is the hardness of getting both span and work low (also known as the transitive-closure bottleneck). Intuitively, for a frontier-based algorithm, the tradeoff embodies in the choice of the frontier size. Choosing a large frontier means to select the vertices for relaxation more aggressively. More vertices can be processed in parallel, but likely much of the work is redundant. We say a vertex $v$ is settled if its tentative distance equals to the real shortest distance, and unsettled otherwise. If an unsettled vertex is used to relax its neighbors, then all work is \emph{ineffective}.
Having more vertices processed in parallel results in more unsettled vertices being processed.
An extreme case of large frontier is the parallel Bellman-Ford algorithm where the frontier is always all vertices. The algorithm will finish in $O(T_{\mathtext{sssp}})$ {step}{s}, where $T_{\mathtext{sssp}}$ is the shortest path tree depth. We note that without shortcutting, this is the lower bound for the number of {step}{s}. However the total work is $O(mn)$.
Especially when the graph is sparse and the shortest path tree is deep, this overhead is significant. For example, in most road networks with millions of vertices, the number of effective updates in each {step}{} is usually only a few hundreds, and thus the extra work is unaffordable.
On the other hand, using a small frontier selects the vertices for relaxation more conservatively, which avoids ineffective work. However it increases the number of {step}{s} needed. Once a vertex is settled, we want to process it as soon as possible, instead of waiting until it is on the top of the priority queue. An extreme case is the parallel Dijkstra's algorithm where the frontier is always just the closest vertex. While no redundant work is done, the algorithm suffers from $\Omega(n)$ span. This long dependence is more significant especially in dense graphs where the shortest path tree is shallow. In fact, in most social networks with billions of vertices, the parallel Bellman-Ford finishes in no more than 50 rounds. Therefore, only using a small frontier delays many effective updates, and harms parallelism.

We plan to tackle this problem, both in theory and in practice. By carefully inspecting and understanding previous algorithms, we propose several simple conditions $\mathcal{I}_{{\mathtext{ex}}}$ and $\mathcal{I}_{{\mathtext{add}}}$ for {Soft Batch-PQ}{}. Some initial experiments show that our algorithm is efficient, and outperforms many existing algorithms (see details below). We plan to study the theory behind the practical results. If successful, we believe we can achieve the first practical parallel SSSP\xspace{} algorithm with non-trivial worst-case cost bounds.
In summary, our goals include the following.

\begin{itemize}
  \item Design the interface of a parallel Batch-PQ\xspace{} and its implementation.
  \item Design a parallel SSSP\xspace{} that uses Batch-PQ\xspace{}, that is practical and theoretically efficient (provable worst-case work and span bounds). We hope that our algorithm could avoid using structure-based parameter. If necessary, parameters could depend on the number of vertices and edges.
  \item Apply the SSSP\xspace{} algorithm to more algorithms. Currently we believe the same data structure applies to s-t shortest path, $A^*$ search, and possibly Huffman tree construction.
  \item For an original algorithm that uses priority queues, formalize the condition that a Batch-PQ\xspace{} data structure could be used to parallelize the original problem. This includes both the correctness guarantee and possibly cost bounds.
\end{itemize}

\textbf{Preliminary Experimental results.} Based on the framework in Algorithm \ref{algo:sssp}, we designed two simple heuristics of the priority queue. The core idea is to keep the simplicity and practicability of the $\Delta$-Stepping\xspace{} algorithm, while absorbing the merits in the theoretically-efficient algorithms.
The two strategies, noted as $\rho$-Stepping\xspace{}-Dense and $\rho$-Stepping\xspace{}-Sparse, are designed for dense and sparse graphs respectively.
The initial frontier size ${\kappa}$ is a fixed size ${\kappa}$ of the current queue.
As mentioned, in a dense version, the number of {step}{s} is likely small. In this case we can afford processing many vertices in a {step}{} to let the tentative distance drop quickly.
To avoid too much redundant work, we also assign closer vertices higher priorities to let them be processed in earlier {substep}{s}.
Therefore in the current $Q$, we identify $d_j$ as the distance of the $j\cdot {\kappa}$-th element smallest. These distance thresholds are used as virtual ``barriers'' to make use of the priority information to roughly order the vertices in the frontier.
In a sparse version, also as mentioned, the main goal is to reduce ineffective work, and thus we focus on first settling all vertices within distance $d_1$, which is some threshold based on the current frontier. We implemented these two simple strategies, each taking about 200 lines of code, to existing SSSP\xspace{} algorithms in highly-optimized graph processing systems, including Ligra \cite{}, Julienne \cite{}, GAPBS, and Galois. We use different graphs types, with up to billions of vertices and edges. Our light-weighted implementation outperforms xx in xxx and xx in xx. We believe our current version is highly practical. We also believe that by further studying the theory, our implementation could get even better performance.

\textbf{Theoretical Study Directions.} As mentioned, theoretically-efficient SSSP\xspace{} algorithm is hard, and reducing both work and depth of the problem is conjectured to be impossible. As mentioned, without adding shortcuts, the lower bound of span is the shortest distance tree depth. However adding shortcuts could be impractical because the number of needed shortcuts is usually much more than the input edges. We propose to study the parameterized bounds. This means that the bound is possibly tied to some properties of the input graph. One example is to bound the span in the shortest path tree depth of the graph. Another formalization we believe useful is the $(k,\rho)$-graph proposed by the PIs in 2016 \cite{}. Informally, a $(k,\rho)$-graph means that any vertex can reach its $\rho$-nearest neighbors in $k$-hops. This basically says that by $k$ {substep}{s}, a settled vertex can settle another $\rho$ vertices. This concept has been used to show the theoretical results in Radius-Stepping\xspace{} \cite{}. We believe it is the key factor that affects how much parallelism we can achieve in an SSSP\xspace{} algorithm. We would like to provide analysis and worst-case bound based on the concept of $(k,\rho)$-graph.

\textbf{Generalization. } If successful, we believe our approach is applicable to other algorithms, including s-t shortest path, $A^*$, and Huffman tree construction. We plan to study the the condition that a {Soft Batch-PQ}{} data structure could be used to parallelize the original problem using a regular \pq{}. This includes both the correctness guarantee and cost bounds.

\begin{algorithm}[t]
\caption{Parallel SSSP algorithm framework using soft priority queue}
\label{algo:sssp}
\SetKwProg{mystruct}{struct}{}{}
\SetKwProg{myfunc}{function}{}{}
\SetKwFor{ParForEach}{parallel foreach}{do}{endfch}
\KwIn{Weighted input graph $G(V, E)$, source $s$}
\KwOut{Shortest distance $d(\cdot)$ of all vertices to $s$}
  \vspace{.5em}
Initialize $d(s)=0$, and all other $d(v)=\infty$ for $v\ne s$. Add $s$ to SoftPQ $\mathtext{\pq}$\\
\While(\tcp*[f]{\textbf{{step}{s}}}) {$\mathtext{\pq}$ is not empty} {
  $\mathcal{F}\gets\mathtext{\pq}$.extract($\mathcal{I}_{{\mathtext{ex}}}$) \tcp*{get the initial frontier}

  \lWhile(\tcp*[f]{\textbf{{substep}{s}}}.) {$\mathcal{F}'\ne \emptyset$} {
    $\mathcal{F} \gets \pq$.relax\_and\_add($\mathcal{F}, \mathcal{I}_{{\mathtext{add}}}$
  }
}
\vspace{.5em}
\mystruct{$\text{SoftPQ}$}{
  Let $Q$ be the set of all elements in this soft priority queue\\
  \myfunc{$\mathcal{F}=$extract($\mathcal{I}_{{\mathtext{ex}}}$)} {
    $\mathcal{F} \gets $All vertices in $Q$ with condition $\mathcal{I}_{{\mathtext{ex}}}$\\
    Remove $\mathcal{F}$ from $Q$\\
  }
  \vspace{.2em}
  \myfunc{$\mathcal{F}'=$relax\_and\_add($\mathcal{F}, \mathcal{I}_{{\mathtext{add}}}$)} {
    $\mathcal{F}' \gets \emptyset$\\
    \ParForEach {$v\in \mathcal{F}$} {
        \ParForEach{$u$ is neighbor of $v$} {
          \If {$v$ relaxes $u$} {
            \lIf {$\mathcal{I}_{{\mathtext{add}}}(u)$} {Add $u$ to $\mathcal{F}'$. Remove $u$ from $Q$ if $u\in Q$}
            \lElse{add $u$ to $Q$ if $u\notin Q$}
          }
        }
    }
  }
}
\end{algorithm}

\hide{
\begin{table}
  \centering
  \small
  \begin{tabular}{llll}

    \bf Algorithm & $\boldsymbol{\mathcal{I}_{{\mathtext{ex}}}(u)}$ & $\boldsymbol{\mathcal{I}_{{\mathtext{add}}}(u)}$ & $\boldsymbol{C_{\mathtext{stop}}}$ \\
    \hline
    Dijkstra's & $d(u)=\min_{v\in Q}$ & \emph{false} & $\mathcal{F}\neq \emptyset$\\
    Bellman-Ford & \emph{true} & \emph{false} & $\mathcal{F}\neq \emptyset$\\
    $\Delta$-Stepping\xspace{} \cite{} & $d(u)\le i\Delta$, $i$ is the current step&$d(u)\le i\Delta$, $i$ is the current round& $\mathcal{F}\neq \emptyset$\\
    Radius-Stepping\xspace{} \cite{} & $d(u)\le r_i$, $r_i$ is the current radius&$d(u)\le r_i$, $r_i$ is the current radius& $\mathcal{F}\neq \emptyset$\\
    Shi and Spencer \cite{}&$d(u)$ is top $k$ in $Q$, where $k$ is a & $d(u)$ is top $k$ in $Q\cup \mathcal{F}'$, and $u$ is & within 4\\
    &fixed parameter&within 4 hops from initial frontier& {substep}{s}\\
    \hline
    $\rho$-Stepping\xspace{}-Dense& $d(u)$ is in the top $k$ in $Q$, where & $d(u)$ is among the top $j\times k$ in $Q$, &within $h$ \\
    &$k$ is a constant fraction of $|Q|$ & $j$ is the current {substep}{} &{substep}{s} \\
    $\rho$-Stepping\xspace{}-Sparse& $d(u)$ is in the top $k$ in $Q$, where& Return \emph{true}{} if $d(u)\le d'$ and $u$ is &within $h$ \\
    &$k=c|Q|$ for constant $c$&  within $h$-hops from the initial&{substep}{s} \\
    &Let $d'=\max_{u}d(u)$ for $u$ in $\mathcal{F}$ & frontier. Return \emph{false}{} otherwise.& \\
    \hline
  \end{tabular}
  \caption{a}\label{tab:allsssp}
\end{table}
}

\hide{
\begin{table}
  \centering
  \small
  \begin{tabular}{lll}

    \bf Algorithm & $\boldsymbol{\mathcal{I}_{{\mathtext{ex}}}(u)}$ & $\boldsymbol{\mathcal{I}_{{\mathtext{add}}}(u)}$ \\
    \hline
    Dijkstra's & $d(u)=\min_{v\in Q}$ & \emph{false}\\
    Bellman-Ford & \emph{true} & \emph{false}\\
    $\Delta$-Stepping\xspace{} \cite{} & $d(u)\le i\Delta$, $i$ is the current step&$d(u)\le i\Delta$, $i$ is the current round\\
    Radius-Stepping\xspace{} \cite{} & $d(u)\le r_i$, $r_i$ is the current radius&$d(u)\le r_i$, $r_i$ is the current radius\\
    Shi and Spencer \cite{}&$d(u)$ is top ${\kappa}$ in $Q$, where $k$ is a & $d(u)$ is top ${\kappa}$ in $Q\cup \mathcal{F}'$ and $j\le 4$, \\
    &fixed parameter&where $j$ is the current round of {substep}{}\\
    \hline
    $\rho$-Stepping\xspace{}-Dense& $d(u)$ is top ${\kappa}_0$ in $Q$, where & $d(u)$ is top $j\times {\kappa}_0$ in $Q$ and $j\le h$, \\
    &${\kappa}={\alpha}|Q|$ for constant ratio ${\alpha}_D$ & where $j$ is the current round of {substep}{}, \\
    &&and $h_D$ is a fixed parameter.\\
    $\rho$-Stepping\xspace{}-Sparse& $d(u)$ is in the top ${\kappa}_S$ in $Q$, where& Return \emph{true}{} if $d(u)\le d'$ and $j\le h_S$, where\\
    &${\kappa}_S={\alpha}_S|Q|$ for constant ratio ${\alpha}_S$& $j$ is the current round of {substep}{}, and $h_S$ is a \\
    &Let $d'=\max_{u}d(u)$ for $u$ in $\mathcal{F}$ &fixed parameter. Return \emph{false}{} otherwise. \\
    \hline
  \end{tabular}
  \caption{a}\label{tab:allsssp}
\end{table}}

\begin{table}
  \centering
  \small
  \begin{tabular}{lll}

    \bf Algorithm & $\boldsymbol{\mathcal{I}_{{\mathtext{ex}}}(u)}$ & $\boldsymbol{\mathcal{I}_{{\mathtext{add}}}(u)}$ \\
    \hline
    Dijkstra's & $d(u)=\min_{v\in Q}$ & \emph{false}\\
    Bellman-Ford & \emph{true} & \emph{false}\\
    $\Delta$-Stepping\xspace{} \cite{} & $d(u)\le i\Delta$, $i$ is the current step&$d(u)\le i\Delta$, $i$ is the current round\\
    Radius-Stepping\xspace{} \cite{} & $d(u)\le r_i$, $r_i$ is the current radius&$d(u)\le r_i$, $r_i$ is the current radius\\
    Shi and Spencer \cite{}&$d(u)$ is top ${\kappa}$ in $Q$, where ${\kappa}$ is a & $d(u)$ is top ${\kappa}$ in $Q\cup \mathcal{F}'$ and $j\le 4$, \\
    &fixed parameter&where $j$ is the current round of {substep}{}\\
    \hline
    $\rho$-Stepping\xspace{}-Dense& $d(u)<d_1$, where $d_j$ is set to be the &
    $u\in \mathcal{F}_0$ and $d(u)<d_j$ and $j\le h$, where $j$ is\\
    &distance of the $(j\cdot {\kappa}_D)$-th smallest  &
     the current round of {substep}{}, $\mathcal{F}_0$ the initial \\
    &element in current $Q$. ${\kappa}_D$ is a constant&
    frontier, and $h_D$ a fixed parameter.\\
    $\rho$-Stepping\xspace{}-Sparse& $d(u)<d_1$, where $d_1$ is set to be the &
    Return \emph{true}{} if $d(u)\le d_1$ and $j\le h_S$, where\\
    &distance of the ${\kappa}_S$-smallest element&
    $j$ is the current round of {substep}{}, and $h_S$ a \\
    & in current $Q$. ${\kappa}_S$ is a constant&
    fixed parameter. Return \emph{false}{} otherwise. \\
    \hline
  \end{tabular}
  \caption{a}\label{tab:allsssp}
\end{table}

\section{Parallel Batch-Dynamic Priority Queue}
\label{sec:sssp}
Another typical example of the abstract data type (ADT) is the priority queue (\pq{}).
A priority queue maintains a set of elements with keys as priorities, and supports three primitives: \mf{Extract-Min} that pops the element with the smallest key, \mf{Insert} that adds a new element, and \mf{Dec-Key} that decreases the key of an element.
Other operations such as \mf{Delete} or \mf{Update} can be implemented via these primitives.
Priority queues have many applications, such as in Prim's algorithm~\cite{prim1957shortest} for computing minimum spanning tree, constructing a Huffman tree~\cite{huffman1952method}, and best-first searches~\cite{pearl1984intelligent}.
The most famous application of priority queues is probably in Dijkstra's algorithm for single-source shortest paths (SSSP\xspace{}) on a graph.
Sequentially, Dijkstra's algorithm uses a priority queue and it is simple and practical.
Similar to the previous section, we will use SSSP and Dijkstra as a motivating example of parallel priority queues.

As opposed to many other graph algorithms, SSSP is notoriously hard in the parallel setting.
To date, there is no work-efficient solution with a non-trivial ($o(n)$) span bound, where $n$ is the number of vertices.
Directly parallelizing Dijkstra's leads to $\Omega(n)$ span due to the use of priority queue: \mf{Extract-Min} pops one vertex at a time for $n$ times in total.
Thus, it is natural to ask if we can parallelize the priority queue.

\myparagraph{Parallelizing priority queue.} There are two major approaches for parallelizing the priority queue.
One is to allow concurrent updates to the priority queue~\cite{alistarh2015spraylist,sundell2005fast,linden2013skiplist,shavit2000skiplist,liu2012lock,henzinger2013quantitative}, and the implementations are based on skip-lists.
The other is to relax the priority queue to allow batch updates and batch extracts on $k$ elements, so multiple processors can work on different vertices simultaneously.
To support these functions, the underlying data structure is a search tree~\cite{blelloch2016parallel}, which can be either the parallel search trees by the PI Sun and coauthors~\cite{Blelloch2016just,sun2018pam,sun2019parallel} in the shared-memory setting, or B-trees~\cite{bingmann2015bulk,sanders2019sequential} in the distributed and external-memory settings.
Unfortunately, maintaining either a skip-list or a search tree is very costly because the total ordering of all elements is kept, which is not required in the sequential heap implementations.

\myparagraph{Practical approaches for SSSP.} Due to the aforementioned problems, practical parallel SSSP algorithms \emph{do not} use the parallel priority queue abstraction.
In practice, Bellman-Ford~\cite{bellman1958routing,ford1956network} can easily be parallelized by relaxing all edges in one round simultaneously, and it performs very well on small-diameter graphs~\cite{ShunB2013,zhang2018graphit}.
Bellman-Ford does not use a priority queue.
Inspired by the good performance of Bellman-Ford, practical parallel SSSP implementations~\cite{zhang2020optimizing,dhulipala2017,nguyen2013lightweight,beamer2015gap} are based on $\Delta$-stepping~\cite{meyer2003delta}, which is a generalization of Bellman-Ford.
In practice, these implementations are faster than those using priority queues since they use flat arrays but not a tree-based priority queue.
However, given a graph with $n$ vertices and $m$ edges, $\Delta$-stepping requires $O(nm)$ work and $O(n)$ span in the worst case, which is no better than the simpler Bellman-Ford.
Also, the parameter $\Delta$ in $\Delta$-stepping highly depends on the input graph structure and the distribution of edge weights, and its best choice is not known without exhaustive searching on each input graph instance~\cite{zhang2020optimizing}.

\myparagraph{Theoretical approaches for SSSP.}
Given the importance of SSSP, there has been a rich literature of theoretical parallel algorithms~\cite{ullman1991high,klein1997randomized,cohen1997using,Shi1999,cohen2000polylog,spencer1997time,blelloch2016parallel} that uses $o(nm)$ work and $o(n)$ span.
Unfortunately, all of these algorithms have some issues that make them potentially impractical on real-world large-scale graphs, including the latest, conceptually simpler algorithm Radius-Stepping\xspace{}~\cite{blelloch2016parallel} proposed by both PIs and coauthors in 2016 (see more details later).
All of these algorithms~\cite{blelloch2016parallel,Shi1999,ullman1991high,klein1997randomized,cohen1997using} relies on some assumptions which make them potentially impractical.
For instance, they all use shortcuts or hopsets to get the corresponding bounds.
To achieve $O(n^{1-\epsilon})$ span in these algorithms, we need to add $\Omega(m^{1+\epsilon})$ shortcut edges.
In practice, the memory size of a single machine can only store the edges in large real-world graphs, and cannot afford to store asymptotically more shortcut edges.
As a result, to date, the problem for exact parallel SSSP remains open:
\begin{center}
  \emph{Q2: Is there a parallel SSSP algorithm that is both practically efficient (faster than Bellman-Ford and $\Delta$-stepping in performance), and theoretically efficient (at least have better bounds than Bellman-Ford)?}
\end{center}

We propose to answer \emph{Q2} by combining existing theoretical and practical results.
On the theory side, the Radius-Stepping\xspace{} paper~\cite{blelloch2016parallel} proposed a concept of $(k,\rho)$-graph, which we believe is useful to analyze parallel SSSP algorithms.
On the practical side, we would like to take implementation details into consideration, and in particular, to avoid using expensive shortcuts or hopsets, or tree-based data structure.

\myparagraph{Proposed work: Parallel Batch-Dynamic Priority Queue (Batch-PQ\xspace{})}.
We believe that \emph{Q2} is a long-standing open problem is due to the lack of a good abstract data type (ADT) of the parallel priority queue. A good abstraction of priority queue should be able to unify the three aspects as discussed above: theory research on parallel SSSP algorithm design and analysis, data structure research for parallel priority queue, and system/software research on parallel SSSP implementation.
Unfortunately, existing work on parallelizing priority queues~\cite{alistarh2015spraylist,sundell2005fast,linden2013skiplist,shavit2000skiplist,liu2012lock,henzinger2013quantitative,bingmann2015bulk,sanders2019sequential} are only used in Dijsktra, not to other theoretical or practical approaches;
existing theoretical approaches~\cite{ullman1991high,klein1997randomized,cohen1997using,Shi1999,cohen2000polylog,spencer1997time,blelloch2016parallel} do not explicitly use parallel priority queues and do not have experimental evaluations;
existing practical implementations~\cite{zhang2020optimizing,dhulipala2017,nguyen2013lightweight,beamer2015gap} are done by the software and system community and they only optimize $\Delta$-stepping, not other algorithms and data structures.
Given the difficulty of parallel SSSP, it seems less likely that a single solution can magically solve all challenges.
We hope our proposed work can help to combine and connect the three aspects in parallel SSSP research.
We believe defining the following abstraction will be very helpful.

\begin{definition}[Batch-PQ]
For a data structure $\mathcal{Q}\xspace{}$ implementing the Batch-PQ\xspace{}, it stores a set of keys for elements in $S$, and supports the following three operations:

\begin{itemize}[topsep=1.5pt, partopsep=0pt,itemsep=0pt]
  \item $\mf{Extract}(\textsc{Cond}\xspace)$: extract the set of $s\in S$ such that $\textsc{Cond}\xspace(s)$ is true for a user-defined predicator $\textsc{Cond}\xspace$,
  \item $\mf{B-Insert}(S')$: batch insertion, which updates $\mathcal{Q}\xspace$ to contain $S\cup S'$,
  \item $\mf{B-Update}(S', \textsc{Relax}\xspace{})$: batch update, which updates the keys of elements in $S'$ using function \textsc{Relax}\xspace{}.
\end{itemize}
\end{definition}

At first glance, it may seem that the Batch-PQ\xspace{} interface is similar to the previous parallel priority queue~\cite{bingmann2015bulk,sanders2019sequential}, which can extract or update a batch of $k$ elements instead of one element.
However, we note that the main difference is that we allow a user-specific predicator $\textsc{Cond}\xspace$ to extract the elements, which supports a much stronger functionality.
For instance, Batch-PQ\xspace{} can be used to implement parallel Bellman-Ford with $\textsc{Cond}\xspace(v)\equiv \textsc{True}\xspace$ and to relax their neighbor vertices using $\mf{B-Insert}$ and $\mf{B-Update}$, but it seems hard with the interface of the previous parallel priority queue (e.g., $k\mf{-Extract}$ and $k\mf{-Update}$).
Although we omit the details of previous algorithms due to space limitations, we note that the interface for Batch-PQ\xspace{} can easily be used to implement practical SSSP algorithms such as $\Delta$-stepping, and theoretical SSSP algorithms such as Radius-Stepping\xspace{}~\cite{blelloch2016parallel} and Shi and Spencer's~\cite{Shi1999}, albeit requiring light modifications.

Therefore, with the Batch-PQ\xspace{} interface, all three previously isolated areas are connected.
The implementation tricks such as the smart use of flat arrays in~\cite{zhang2020optimizing,dhulipala2017,nguyen2013lightweight,beamer2015gap} can inspire the implementation for theoretical algorithms since they share the same underlying implementation for Batch-PQ\xspace{}.
The theoretical analysis in~\cite{blelloch2016parallel,Shi1999} still holds once we plug in the cost bounds for the Batch-PQ\xspace{} operations.
Also, when we design or implement parallel SSSP algorithms, we can treat the Batch-PQ\xspace{} as a black box, which greatly simplifies the thought process; however, that was not the case when designing previous algorithms either theoretically~\cite{ullman1991high,klein1997randomized,cohen1997using,Shi1999,cohen2000polylog,spencer1997time,blelloch2016parallel} or practically~\cite{zhang2020optimizing,dhulipala2017,nguyen2013lightweight,beamer2015gap}.
We believe the Batch-PQ\xspace{} interface bridges together the theoretical algorithms and analysis, data structure design, and implementation details.
Therefore, we are optimistic about answering \emph{Q2} and designing parallel SSSP algorithms that are efficient both theoretically and practically.

\textbf{Theoretical Study Directions.}
Aiming at bridging the gap between theory and practice, we combined some theory insights in Radius-Stepping\xspace{} and some ideas in practical implementations of $\Delta$-Stepping\xspace{}~\cite{zhang2020optimizing,dhulipala2017}, as well as the Batch-PQ\xspace{} abstraction as proposed above. Bringing them together, we designed a preliminary parallel SSSP algorithm $\rho$-Stepping\xspace{}.
The reason we pick Radius-Stepping\xspace and $\Delta$-Stepping\xspace is that they all run in steps, and in each step, the algorithms ``settle down'' the SSSP distances for vertices in a ``range''.
The difference between Radius-Stepping\xspace and $\Delta$-Stepping\xspace is that, $\Delta$-Stepping\xspace uses a fixed range, while Radius-Stepping\xspace computes a distance range in each step that provides efficient cost bounds.
$\rho$-Stepping\xspace{} adopts the idea similar to Radius-Stepping\xspace that dynamically decides the distance range, and the implementation of Batch-PQ\xspace{} is a flat array-structure similar to $\Delta$-Stepping\xspace.
Such implementation avoids using the tree-base structure and does not add shortcut edges.
As such, $\rho$-Stepping\xspace{} is already very efficient in practice (Figure~\ref{fig:prelim}(b)), but again does not have interesting bounds.

The first question we would like to study is an efficient design for Batch-PQ\xspace.
Our goal is to achieve better bounds than using the flat array in $\Delta$-Stepping\xspace, and the search-tree in Radius-Stepping\xspace that is not work-efficient and practically inefficient due to the pointer-based structure.
Our goal is a hierarchical structure that provides improved bounds, and can be implemented based on arrays that is fast in practice.

The second question is about the bounds for the SSSP algorithm.
Since we cannot add shortcuts in practice, showing interesting worst-case bounds seems much more challenging---there is no much we can do for a chain without adding shortcuts.
However, we believe we can achieve parameterized bounds.
This means that the bounds are possibly tied to some properties of the input graph.
We believe that the concept of $(k,\rho)$-graph proposed by the PIs in 2016~\cite{blelloch2016parallel} can be helpful.
Informally, a $(k,\rho)$-graph means that any vertex can reach its $\rho$-nearest neighbors in $k$-hops.
This basically means that by $k$ {substep}{s}, a settled vertex can further settle $\rho$ vertices.
This concept was used to prove the bounds of Radius-Stepping\xspace{}, and seems to be the key factor on how much parallelism we can achieve in SSSP\xspace{}.
Hence, the PIs propose to explore what parameterized bounds we can show for practical SSSP algorithms, with the tools of Batch-PQ\xspace and $(k,\rho)$-graph.
It seems that showing the bounds for $\rho$-Stepping\xspace is promising due to the similarity to Radius-Stepping\xspace.

\textbf{Preliminary Experimental Results.}
As a proof-of-concept, we implemented $\rho$-Stepping\xspace{} and compared it to state-of-the-art SSSP implementations (see Figure~\ref{fig:prelim}(b)).
We can see that even our preliminary version is already faster than the state-of-the-art graph processing systems on almost all graphs.
This indicates the practicability of our solution.
We will keep exploring better implementations, especially by exploring the parameter space of $\rho$-Stepping\xspace, which seems promising using the Batch-PQ\xspace{} abstraction.
We will also improve the algorithm according to the progress on the theoretical side of the algorithm.

\myparagraph{Other applications using Batch-PQ\xspace{}.}
We note that many problems with sequential algorithms using priority queues remain unsolved in the parallel setting.
We believe the reason is the lack of an elegant interface of the parallel priority queue and efficient underlying implementation.
If the previously proposed work is successful, we can consider Batch-PQ\xspace as a black box and apply it to many other algorithms.

One possible topic to study is other shortest-paths problems such as $s$-$t$ shortest path and $A^*$ search.
There are many textbook algorithms to solve these problems sequentially~\cite{CLRS} by using priority queues.
If the design of Batch-PQ\xspace{} is successful, we want to explore whether we can solve these problems in parallel.

Another interesting question to explore is when the parallel priority queue can be used.
For instance, we cannot directly use Batch-PQ\xspace{} to parallelize Prim's algorithm since we can only guarantee the smallest edge in the priority queue is in the MST, and not even for the second smallest.
The reason we can use Batch-PQ\xspace{} in SSSP is that SSSP allows ``priority-inversion''~\cite{alistarh2019efficiency}.
However, we believe this property is too strong.
In fact, we believe Batch-PQ\xspace{} can be used in more applications such as constructing Huffman tree, which is another typical use case of priority queue but does not allow priority-inversion.
We plan to study the theoretical question about the accurate boundary when Batch-PQ\xspace{} can apply correctly.
\subsection{Comparisons and Discussions}\label{sec:discussion}

For $\rho$-Stepping\xspace, the number of total steps is $O(k_\rho n/\rho)$ for undirected graphs and $O(k_n n/\rho)$ for directed graphs.
The undirected case is a factor of $O(\log \rho L)$ better than Radius-Stepping\xspace (Radius-Stepping\xspace does not have non-trivial span bound on direct graphs).
The work bound is off by a factor of $k_n/k_\rho$ on undirected graphs, but it applies to directed graphs.
Also, our experiments show that, on social and web graphs, $k_n/k_\rho$ is usually small (\cref{fig:k-rho}) for reasonably large values of $\rho$ (e.g., $\rho > \sqrt{n}$).

Both $\rho$-Stepping\xspace and $\Delta^*$-Stepping\xspace focus on practical considerations.
Since in practice we usually pick a large $\rho$, the number of steps is small. This leads to a small overhead for step-based synchronization.
\cref{thm:deltass} show that $\Delta^*$-Stepping\xspace only incurs a factor of $1+L/\Delta$ more steps (recall $L=\max w(e)$) than Bellman-Ford, upper bounding the synchronization cost in practice (\cref{fig:visted-per-step}).
Regarding work, \cref{thm:rho-steps,thm:deltass} show that both tournament tree\xspace{}-based and array-based versions are efficient when using proper parameters of $\rho$ and $\Delta$.
Exactly in our experiments, the best values $\rho$ and $\Delta$ match the analysis here (e.g., a large $\rho$ on social networks).
We note that Bellman-Ford has better work and span than both $\rho$-Stepping\xspace and $\Delta^*$-Stepping\xspace.
In fact, it seems hard to beat the work and span of Bellman-Ford (parameterized on $k_n$) if no shortcut edges are allowed.
Our analysis provides worst-case guarantees for $\rho$-Stepping\xspace and $\Delta^*$-Stepping\xspace, and they seem good for the $(k_{\rho},\rho)$ parameters of many real-world graphs.
In practice, both $\rho$-Stepping\xspace{} and $\Delta^*$-Stepping\xspace{} exhibit better performance than Bellman-Ford because of visiting fewer vertices and edges (more efficient ``work'').
Since analyzing SSSP algorithms based on $(k,\rho)$-graph\xspace is new, many interesting questions remain open.

The work for Radius-Stepping\xspace and Shi-Spencer can be improved by at most a logarithmic term.

\section{Experiments}
\label{sec:exp}

\newcommand{scale-free network}{scale-free network}
\input{tables/overalltable.tex}

\myparagraph{Experimental setup.}
\label{sec:exp:setup}
We run all experiments on a quad-socket machine with Intel Xeon Gold 6252 CPUs with a total of 96 cores (192 hyperthreads). The system has 1.5TB of main memory and 36MB L3 cache on each socket. Our codes were compiled with \texttt{g++} 7.5.0 using \texttt{CilkPlus} with \texttt{-O3} flag. For all parallel implementations, we use all cores  and \texttt{numactl -i all}, which evenly spreads the memory pages across the
processors in a round-robin fashion.

We implemented three algorithms based on the framework in \cref{sec:framework}: Bellman-Ford (\textit{PQ-BF}\xspace), $\Delta^*$-Stepping\xspace{} (\textit{PQ-$\Delta^{*}$}\xspace), and $\rho$-Stepping\xspace{} (\textit{PQ-$\rho$}\xspace).
We use array-based {\textsc{LaB-PQ}}\xspace{} because we observe that when the output size of \mf{Extract}\xspace{} is large, the array-based implementation has better performance than tournament tree\xspace{} (see more details in the full paper).
For all graphs we use, the best running time is achieved using a reasonably large $\rho$.
We compare our implementations with state-of-the-art SSSP implementations: Bellman-Ford algorithm in Ligra~\cite{ShunB2013}, $\Delta$-Stepping\xspace{} in Julienne~\cite{dhulipala2017}, GAPBS~\cite{beamer2015gap,zhang2020optimizing}, and Galois~\cite{nguyen2013lightweight}. Throughout the section, when we refer to ``$\Delta$-Stepping\xspace{}'', it includes our $\Delta^*$-Stepping\xspace{}, and the existing $\Delta$-Stepping\xspace{} in Julienne, GAPBS and Galois.

We test seven graphs, including four social networks com-orkut (OK)~\cite{yang2015defining}, Live-Journal (LJ)~\cite{backstrom2006group}, Twitter (TW)~\cite{kwak2010twitter} and Friendster (FT)~\cite{yang2015defining}, one web graph WebGraph (WB) \cite{webgraph}, and two road graphs~\cite{roadgraph} RoadUSA (USA) and Germany (GE). The graph information is provided in Table \ref{tab:alltime}.
In almost all experiments, the social and web graphs show a similar trend.
This is because they follow similar power-law-like degree distribution. Throughout the section, we use ``scale-free network{s}'' to refer to social and web graphs.
On scale-free network{s}, we set edge weight uniformly at random in range $[1,2^{18})$. On road graphs, the edge weights are from the original dataset, which is up to $2^{25}$.

For all $\Delta$-Stepping\xspace{} algorithms (except for \cref{fig:intro:delta} where we vary~$\Delta$), we report the \emph{best running time} across all $\Delta$ values. When we report average of multiple sources, we first find the best $\Delta$ value on one source, and use it for other sources. We do this for every graph-implementation combination.
For all $\rho$-Stepping\xspace{} algorithms (except for in \cref{tab:alltime} where we explicitly report the best running time across $\rho$ values), we use a fixed value of $\rho=2^{21}$. For most experiments, we report the average of 10 sources.
When taking the average is meaningless, we use one representative source.

In this section, we will first discuss the overall performance of all implementations. We then compare some statistics to better understand the performance of \textit{PQ-$\rho$}\xspace{}, \textit{PQ-$\Delta^{*}$}\xspace{} and \textit{PQ-BF}\xspace{}. We evaluate the number of vertices visited by the algorithm as an indicator of the overall work. Since road graphs exhibit different properties from the scale-free network{s}, we then discuss road graphs separately. We also analyze the two algorithms $\rho$-Stepping\xspace{} and $\Delta$-Stepping\xspace{} with their corresponding parameters.
Due to page limit, we will discuss the $k_\rho$ properties for each graph in the full version, and only show the $k_\rho$-$\rho$ curves in this paper (\cref{fig:k-rho}).
\ifx0\undefined
Due to page limit, we postpone some figures, discussions and more experiments to the full paper~\cite{ssspfull} (e.g., using different machines and different types of sources).
\else \fi
Generally, we show that our $\rho$-Stepping\xspace{} has \textbf{especially good performance on scale-free network{s}}, and the performance gain of $\rho$-Stepping\xspace is from three aspects: \textbf{good parallelism}, \textbf{less overall work}, and \textbf{more evenly distributed work to all steps}. We summarize conclusions and interesting findings at the end of this section.

\myparagraph{Overall Performance.} We present the running time of all
implementations in \cref{tab:alltime}.
In all cases, one of our implementations achieves the best performance, and is 1.14$\times$ to orders of magnitude faster than the previous implementations. We show a heat map of relative parallel running time in \cref{fig:exp:heat}.

On scale-free network{s}, \textit{PQ-$\rho$}\xspace{} and \textit{PQ-$\Delta^{*}$}\xspace{} outperform all existing implementations.
\textit{PQ-$\rho$}\xspace{} has better performance.
On average over five graphs, \textit{PQ-$\rho$}\xspace{} is $1.41\times$ faster than Galois, $1.83\times$ faster than Julienne and GAPBS, and $1.93\times$ faster than Ligra.

On road graphs, \textit{PQ-$\Delta^{*}$}\xspace{} is the fastest, and \textit{PQ-$\rho$}\xspace{} is also competitive.
Ligra did not finish in 30 seconds on road graphs, since Ligra uses plain Bellman-Ford that is inefficient for graph with deep shortest-path tree (more than $10^4$, see \cref{fig:k-rho}).
Our \textit{PQ-BF}\xspace{} with the neighbor-set optimization (see \cref{sec:implementation}) finishes on both graphs in about 0.4s.

We report the sequential running time of the corresponding parallel version and show self-speedup in Table \ref{tab:alltime}.
We note that comparing the sequential running time
of different implementations does not seem useful because both $\Delta$-Stepping\xspace and $\rho$-Stepping\xspace are parameterized.
To get the best sequential performance, one should just use a small $\Delta$ or $\rho$.
The reported time is the sequential performance using the corresponding parameter that performs best in parallel, and it makes more sense just to compare the speedup numbers.
The self-speedup of \textit{PQ-$\rho$}\xspace{} is almost always the best among all implementations (\textit{PQ-$\Delta^{*}$}\xspace{} is close but slightly worse). Hence, the good performance of \textit{PQ-$\rho$}\xspace{}, especially on scale-free network{s}, is partially due to good scalability.  In other words, \textit{PQ-$\rho$}\xspace{} achieves the best ``work-span tradeoff'' in practice.

Among the implementations of the same algorithm, \textit{PQ-BF}\xspace{} outperforms Ligra on all graphs.
For all $\Delta$-Stepping\xspace{} algorithms, \textit{PQ-$\Delta^{*}$}\xspace{} is also the fastest on all graphs.
Overall, our three algorithms outperform existing implementations, indicating the efficiency of stepping algorithm framework for parallel SSSP implementations.

\label{sec:exp:analysis}

\hide{
\begin{figure}[t]
  \centering
  \includegraphics[width=\columnwidth]{figures/heat.pdf}
  \vspace{-.2in}
  \caption{ \small \label{fig:exp:heat} \textbf{Parallel running time relative to the fastest on each graph.} \mdseries Implementations noted with $*$ are ours. ``Ave.'' means average on scale-free network{s} (over the five graphs) and road graphs (over the two graphs), respectively. (96h) means parallel running time (96 cores, 192 hyperthreads). (1) means the sequential running time of the corresponding parallel version. (SU) means the speedup number.
  \textit{PQ-$\rho$}\xspace{}-fix denotes uses a fixed parameter $\rho$ ($2^{21}$ for scale-free network{s} and $2^{15}$ for road graphs). \textit{PQ-$\rho$}\xspace{}-best is the best running time using all values of $\rho$.

  }
  \vspace{-.2in}

\end{figure}
}

\input{expfigures.tex}

\myparagraph{Number of visits to vertices.}
Unlike Dijkstra, other parallel SSSP algorithms can visit each vertex or edge more than once.
While this allows for parallelism, the total work is also increased.
To show how much ``redundant'' work is done for the stepping algorithms, we measure the average number of visits per vertex (\cref{fig:enqueue}), and the number of visited vertices in each step on four representative graphs (\cref{fig:visted-per-step})
\ifx0\undefined
\footnote{We also measured the number of visited edges, which show very similar trend to the vertices. Due to page limit, we report the numbers in the full version of the paper, which also shows the results for all seven graphs.}.
\else
\footnote{We also measured the number of visited edges, which show very similar trend to the vertices. We present the results in \cref{fig:visted-per-step}}.
We show results for four representative graphs, and full results are in \cref{fig:visted-per-step}.
\fi
We note that the other systems vary a lot in implementation details, and it is hard to directly measure these quantities from their code.
Hence, we compare among our implementations.
For the same reason, \textit{PQ-$\Delta^{*}$}\xspace{} may not precisely reflect the numbers of other $\Delta$-Stepping\xspace{} implementations. In this paragraph, we first focus on the scale-free network{s}, and discuss road graphs later.

\Cref{fig:enqueue} shows the average number of visits per vertex.
{{\ifx0\undefined}}
On the two small graphs (OK and LJ), since the work cannot saturate all 192 threads, \textit{PQ-$\rho$}\xspace act similar to Bellman-Ford to maximize parallelism and uses visits more vertices than $\Delta$-Stepping\xspace{}.
\else
On the two small graphs (OK and LJ), since the work cannot saturate all 192 threads, both \textit{PQ-$\Delta^{*}$}\xspace and \textit{PQ-$\rho$}\xspace act similar to Bellman-Ford to maximize parallelism.
\fi
{{\ifx0\undefined}}
For the larger graphs (TW, FT, and WB), \textit{PQ-$\rho$}\xspace{} always triggers the smallest average visit to vertices.
\else
For the three larger graphs (TW, FT, and WB), \textit{PQ-$\rho$}\xspace{} always triggers the smallest average visit to vertices and edges.
\fi
The trend showed in \cref{fig:enqueue} exactly matches the sequential time of each implementation.
Hence, one advantage of \textit{PQ-$\rho$}\xspace{} over \textit{PQ-BF}\xspace{} and \textit{PQ-$\Delta^{*}$}\xspace{} on scale-free network{s} is less total work.

\Cref{fig:visted-per-step} shows the number of visited vertices per step.
In \textit{PQ-BF}\xspace{}, the numbers always grow quickly to a large value, stay for a few steps, and finish quickly.
Although usually using the fewest steps, \textit{PQ-BF}\xspace{} is the slowest, since the dense steps cause many redundant relaxations.
\textit{PQ-$\Delta^{*}$}\xspace{} usually uses more steps than both \textit{PQ-BF}\xspace and \textit{PQ-$\rho$}\xspace{}.
In most of the steps (at the beginning and the end), \textit{PQ-$\Delta^{*}$}\xspace{} visits only a small number of vertices, but the peak values are much higher than \textit{PQ-$\rho$}\xspace{}.
\textit{PQ-$\rho$}\xspace{} shows a more even pattern across the steps: in most of the steps, it processes a moderate number of vertices, and the peak value is much smaller than \textit{PQ-$\Delta^{*}$}\xspace{} or \textit{PQ-BF}\xspace{}.

These patterns in \cref{fig:visted-per-step} reflect the nature of the three algorithms.
Bellman-Ford always visits all vertices in the frontier in each step.
This created significant redundant work.
$\Delta^*$-Stepping\xspace{} controls work-span tradeoff based on distances.
On scale-free network{s}, it reaches the peak work in some middle steps, which is significantly higher than other steps.
$\rho$-Stepping\xspace{} controls the work-span tradeoff using the number of vertices processed per step.
We believe on scale-free network{s}, this quantity is a closer indicator to the actual ``work'' in each step than the distance gap is.
In other words, $\rho$-Stepping\xspace{} controls the work in each step that is minimal to saturate all processors, so it explores sufficient parallelism with minimized redundant work.

\ifx0\undefined
\else
\myparagraph{$k$-$\rho$ property for different graphs.} To better understand the properties of different graphs, we show the $\rho$-$k_\rho$ curves in \cref{fig:k-rho}. We note that computing the exact value of $k_\rho$ is expensive (which makes Radius-Stepping\xspace{} impractical), and we estimate the value using 100 samples.
As mentioned, this indicates how much potential parallelism we can get on the original graph (i.e., without shortcuts). We show the value of $k_\rho$ when $\rho$ is $\log n$, $\sqrt{n}$, $n/\log n$, $n/10$, and $n$, respectively. Note that $k_n$ is the shortest path tree depth.

On all scale-free network{s}, we observe that $k_n$ is $O(\log n)$ (about $2\log n$). Interestingly, all scale-free network{s} are $(\log n,\sqrt{n})$-graph, which means that almost all vertices can reach their $\sqrt{n}$ nearest vertices by $\log n$ hops.
This is not surprising for scale-free network{s}, which have some ``hubs'' that are well connected to other vertices.
These vertices are easy to be reached by any source in a few steps.
Once reached, they can reach a lot of other vertices, which quickly accumulate $\sqrt{n}$ nearest neighbors for any source.
The road graphs have a different $\rho$-$k_\rho$ pattern.
On GE and USA, it takes more than $100$ hops to reach $\sqrt{n}$ nearest vertices, and $O(\sqrt{n})$ steps to reach $O(n)$ nearest vertices.
We believe the result is reasonable since road graphs are (almost) planar graphs.
\fi

\myparagraph{Discussions for Road Graphs.}
Road graphs are (almost) planar and have different $k$-$\rho$ patterns than other graphs, and the shortest-path trees are deep and slim.
Hence, without the special optimizations (e.g., in Ligra and Julienne), the performance is slow.
As mentioned in \cref{sec:implementation}, our optimization expands multiple levels in the shortest-path tree in one step.
This makes the performance of our implementations competitive or better than GAPBS and Galois.

On road graphs, \textit{PQ-$\Delta^{*}$}\xspace{} is the fastest.
This somehow indicates that expanding with distance may be a good strategy for road graphs.
One possible reason is that they are planar graphs with Euclidean distance.
Hence, setting fixed-width ``annuli'' seems a reasonable work-parallelism tradeoff, when using a proper $\Delta$.
Since the frontier on road graphs is small, \textit{PQ-$\rho$}\xspace{} has insufficient frontier size in each step for enough parallelism.
Hence, it is hard for \textit{PQ-$\rho$}\xspace{} to control the number of vertices visited precisely, and the performance is slightly slower than \textit{PQ-$\Delta^{*}$}\xspace.
However, \textit{PQ-$\rho$}\xspace{} has more stable performance than \textit{PQ-$\Delta^{*}$}\xspace{} in the parameter space (\cref{fig:intro:delta,fig:rho-time}).

\myparagraph{$\Delta$-Stepping\xspace and $\Delta$.} We test all $\Delta$-Stepping\xspace{} algorithms with varying $\Delta$ on all graphs.
For each test case, we normalize the running time to the best time across all $\Delta$ values.
For page limit, we present four graphs in \cref{fig:intro:delta}, and the full results in the full paper~\cite{ssspfull}.

On the same graph, the best choice of $\Delta$ varies a lot
for different systems.
On TW, Julienne's best $\Delta$ is $2^{12}\times$ larger than Galois's.
The best $\Delta$ for one system can make another system up to $4\times$ slower.
The selection of $\Delta$ in one system does not generalize to other systems.
Secondly, even though all scale-free network{s} have the same edge weight distribution, for the same implementation, the best choice of $\Delta$ varies a lot on different graphs.
Therefore, the selection of $\Delta$ on one graph does not generalize to other graphs.
On the same graph, the performance is sensitive to the value of $\Delta$.
Usually, $2$--$4\times$ off may lead to a 20\% slowdown, and $4$--$8\times$ off may lead to a 50\% slowdown.
A badly-chosen $\Delta$ can largely affect the performance.
As a result, for every graph-implementation combination, we have to search the best parameter $\Delta$. Fortunately, we find out that different sources show relatively stable performance for the same implementation-graph pair. This is also the conventional way of tuning $\Delta$ (we also did so).
\ifx0\undefined
We present the results in the full paper \cite{ssspfull}.
\else
We present the results in \cref{fig:differentsource}.
\fi

{{\ifx0\undefined}}
\else
Also, for almost all systems, the performance on road graphs is more sensitive to the value of $\Delta$ than on scale-free network{s}.
\fi

\myparagraph{$\rho$-Stepping\xspace{} and $\rho$.} We test $\rho$-Stepping\xspace{} with varying $\rho$.
When $\rho$ is small, the running time increases significantly.
This is also due to the lack of parallelism (similar to when $\Delta$-Stepping\xspace{} uses small $\Delta$).
When $\rho$ gets large, the performance drops by no more than 20\%.
The best choices of $\rho$ are very consistent on different graphs. This is because the choice of $\rho$ in practice depends on the right level of parallelism we want to achieve, instead of the graph structure or edge weight distribution.
As discussed, $\rho$-Stepping\xspace{} distributes work more evenly to each step.
The goal of setting $\rho$ is to enable enough work to exploit full parallelism in each step, but without introducing more redundant work.
The performance on road graphs are less sensitive, probably because the frontier size seldom reaches $\rho$ in road graphs.
We also tested $\rho$-Stepping\xspace{} on various machines. We observe that the best choice of $\rho$ is still relatively consistent among different settings. We will present more results in the full paper.

Generally speaking, using large $\Delta$ or $\rho$ gives better (and more stable) performance than small $\Delta$ or $\rho$ values.
This is not surprising because when these parameters are large, $\Delta$-Stepping\xspace{} and $\rho$-Stepping\xspace{} degenerate to Bellman-Ford that still has reasonable performance on social networks.
When the parameters are small, $\Delta$-Stepping\xspace{} and $\rho$-Stepping\xspace{} both degenerate to Dijkstra and loses parallelism.

\myparagraph{Summary. }
\label{sec:exp-sum}
In summary, our \textit{PQ-$\rho$}\xspace{} generally achieves the best performance on the five scale-free network{s}. On average of the five graphs, \textit{PQ-$\rho$}\xspace{} is 1.41-1.93$\times$ faster better than existing systems.
On the two road graphs, \textit{PQ-$\Delta^{*}$}\xspace{} always has the best performance, which is at least 14\% better than existing systems.
The good performance of \textit{PQ-$\rho$}\xspace{} on scale-free network{s} comes from three aspects. The first is scalability, indicated by the good self-speedup.
Secondly, it visits fewer vertices and edges on large scale-free network{s}, which indicates less overall work.
Lastly, the work is more evenly distributed to each step, such that each step can exploit sufficient parallelism, and also avoid performing ``ineffective'' work to relax the neighbors of unsettled vertices.
This also indicates that on scale-free network{s} with uniformly distributed edge weights, controlling the number of vertices visited per step is a good strategy.
On road graphs with Euclidean distance, $\Delta$-Stepping\xspace{} shows better performance.

Our \textit{PQ-$\rho$}\xspace{} implementation generally shows stable performance across $\rho$ values on all tested graphs.
A fixed $\rho$ almost always gives performance within 5\% off the performance with the best $\rho$.

Finally, on all tested graphs, \textit{PQ-$\rho$}\xspace{} and \textit{PQ-$\Delta^{*}$}\xspace{} are faster than all existing SSSP implementations (except for RoadUSA, \textit{PQ-$\rho$}\xspace{} is 0.01s slower than Galois).
\textit{PQ-BF}\xspace{} is faster than Ligra on all graphs.
This indicates the efficiency of the stepping algorithm framework on implementing and optimizing parallel SSSP algorithms.

\section{Frameworks}\label{sec:framework}

\subsection{The {\textsc{LaB-PQ}}\xspace{} Abstraction}
An abstract data type \emph{\textbf{La}zy-\textbf{B}atched \textbf{P}riority \textbf{Q}ueue\xspace{}}, or \emp{{\textsc{LaB-PQ}}\xspace{}}, denoted as $\mathbb{PQ}\xspace$, maintains a universe of records\xspace{} $({\mathit{id}}\xspace{},k)$, where ${\mathit{id}}\xspace{}\in I$ is the unique \emp{identifier} for this record\xspace{} and $k\in K$ is the \emp{key}.
In some applications, each record\xspace{} also has a \emp{value} $v\in V$.
In this paper, if not specified, we assume an empty value type for simplicity.
{{\ifx0\undefined}}
\else
In many applications, the identifier type $I$ is the natural number set $\mathbb{N}$.
\fi
In all SSSP algorithms in this paper, the identifiers are vertex labels from 1 to $n$.
The total ordering of all keys is determined by a comparison function $<_K:K\times K\mapsto \mathit{Bool}\xspace{}$.
A {\textsc{LaB-PQ}}\xspace{} $Q\in \mathbb{PQ}\xspace$ is associated with a \emp{mapping function\xspace{}} $\delta_{Q}:I\mapsto K$, which maps an identifier to its corresponding key (or key-value) that can change dynamically over time.
With clear context, we omit the subscription $Q$, and use $\delta[{\mathit{id}}\xspace]$ to denote the mapping from ${\mathit{id}}\xspace$ to key.
In the SSSP algorithms of this paper, this mapping function\xspace{} maps each vertex label to its (tentative) distance.
In our implementation, this mapping function\xspace{} is passed to {\textsc{LaB-PQ}}\xspace{} by a pointer to the tentative distance array. More formally, a {\textsc{LaB-PQ}}\xspace{} $\mathbb{PQ}\xspace$ is parameterized on the following:
\vspace{.05in}

\begin{center}
\noindent {\small\begin{tabular}{@{}ll@{}}

  \hline
  $\boldsymbol{I}$ & Unique identifier type \\
  $\boldsymbol{K}$ & Key type \\
  $\boldsymbol{V}$ & (Optional) Value type \\
  $\boldsymbol{<_K}:K\times K\mapsto \mathit{Bool}\xspace$ &  Comparison function on $K$\\
  $\boldsymbol{\delta_Q}:I\mapsto K\times V$ & A mapping from an id to its key (or key-value)\\

  \midrule
\end{tabular}}
\end{center}

A {\textsc{LaB-PQ}}\xspace{} maintains a subset of identifiers in the universe.
It can extract records\xspace{} with (relatively) small keys in parallel based on $\delta[\cdot]$.
The interface of the {\textsc{LaB-PQ}}\xspace{} includes two functions: \mf{Update}\xspace{} and \mf{Extract}\xspace{} (see Table \ref{tab:interface})
We note that these two functions are sufficient for SSSP application. We discuss more functionalities of {\textsc{LaB-PQ}}\xspace{} in
\ifx0\undefined
the full version of this paper~\cite{ssspfull}.
\else
\cref{sec:fully-dynamic}.
\fi

\textbf{\textsc{Update}}{}$({\mathit{id}}\xspace)$ function commits an update to $Q$ regarding the record with identifier ${\mathit{id}}\xspace$.
  It ``notifies'' $Q$ that the new key for this record\xspace{} is now in $\delta[{\mathit{id}}\xspace]$.
  If ${\mathit{id}}\xspace$ is not in $Q$ yet, \mf{Update}\xspace{} inserts it to $Q$.
  Multiple \mf{Update}\xspace{} can be executed concurrently.
We note that the change of the record\xspace{} is embodied in the change of $\delta[{\mathit{id}}\xspace]$, and thus the data structure only needs to know the record\xspace{}'s ${\mathit{id}}\xspace$ to address the modification.
An important observation is that, we do not have to modify $Q$ immediately, but can execute them \emph{lazily}.
These changes make no difference to any other operations on $Q$ before the next \mf{Extract}\xspace{}.
Compared to the classic ``batch-dynamic'' setting, our interface avoids explicitly generating the batch, which simplifies the algorithm and improves performance.

\textbf{\textsc{Extract}}\xspace{}$(\theta)$ returns all identifiers in $Q$ with key $\le \theta$, and then deletes them from $Q$.
Note that the result of \mf{Extract}\xspace{} reflects all previous \emph{modifications} to $Q$, including \mf{Update}\xspace{} functions and deletions from the previous \mf{Extract}\xspace{}.
It then extracts the corresponding record\xspace{s} based on the latest
view of $Q$.
An \mf{Extract}\xspace{} function \emph{cannot} be executed concurrently with other functions (\mf{Update}\xspace{} or another \mf{Extract}\xspace{}).
This is required for {\textsc{LaB-PQ}}\xspace{}'s correctness.
\hide{
Generating the sequence as output does not increase the asymptotical cost of the algorithm, but can be expensive in practice because it increases memory footprint. In our implementations, we use some optimizations to avoid physically generating the sequence (more details in Section \ref{sec:implementation}).}

\myparagraph{Augmenting {\textsc{LaB-PQ}}\xspace{}.} In some applications, we need a ``sum'' (the \emph{augmented value} of type $A$) of all record\xspace{s} (keys and possible values) in the {\textsc{LaB-PQ}}\xspace{} .
We refer to this as $Q.\mf{Reduce}\xspace()$.
This function first map each record in $Q$ to a value of type $A$, and use a binary commutative
and associative operator~$\oplus$ ($(A,\oplus)$ is a commutative monoid) to compute abstract sum of all record\xspace{s} in $Q$ using $\oplus$.

\begin{table}

  \small
  \begin{tabular}{@{}p{2.5cm}@{  }@{  }p{5.5cm}@{}}
    \hline

    $Q.\mf{Update}\xspace({\mathit{id}}\xspace)$: \flushright$\mathbb{PQ}\xspace \times I \mapsto \Box$ &

    Modify the record\xspace{} with identifier ${\mathit{id}}\xspace$ in $Q$ to $\delta[{\mathit{id}}\xspace]$. If ${\mathit{id}}\xspace\notin Q$, first add it to $Q$.\\
    $s=Q.$\mf{Extract}\xspace{}$(\theta)$:\flushright $\mathbb{PQ}\xspace\times K\mapsto \mathit{seq}\xspace$& Return identifiers in $Q$ with keys no more than $\theta$ and delete them from $Q$.\\
    \hline
  \end{tabular}\vspace{-.2in}
  \caption{\small \textbf{Interface of {\textsc{LaB-PQ}}\xspace{}.} \label{tab:interface}}\vspace{-.1in}
\end{table}

\subsection{The Stepping-Algorithm Framework}\label{sec:step-algo}

\input{table-framework}
{
\begin{algorithm}
\caption{ The Stepping Algorithm Framework.\label{alg:framework}}
\fontsize{9pt}{10pt}\selectfont
\SetKwFor{parForEach}{ParallelForEach}{do}{endfor}
\KwIn{A graph $G=(V,E,w)$ and a source node $s$.}
\KwOut{The graph distances $d(\cdot)$ from $s$.}

\DontPrintSemicolon
$\delta[\cdot]\leftarrow +\infty$, associate $\delta$ to a {\textsc{LaB-PQ}}\xspace{} $Q$\\
$\delta[s]\leftarrow 0$, $Q.\mf{Update}\xspace{}(s)$\\
\While{$|Q|>0$\label{line:while-loop}} {
\parForEach {$u\in Q.\mf{Extract}\xspace(\mf{ExtDist}\xspace)$\label{line:extract}} {
    \parForEach{$v\in N(u)$\label{line:innerloop2}} {
    \If {$\mf{WriteMin}\xspace(\delta[v],\delta[u]+w(u,v))$ \label{line:writemin}} {
      $Q.\mf{Update}\xspace{}(v)$
    }
  }
}
Execute \finishcon
}
\Return {$\delta[\cdot]$}\\[.15em]
\end{algorithm}
}

On top of the {\textsc{LaB-PQ}}\xspace interface, we also propose a simple \emp{stepping algorithm framework},
in order to reveal the internal connection of the existing SSSP algorithms.
Recall the two sequential textbook algorithms, Dijkstra's algorithm~\cite{dijkstra1959} and Bellman-Ford algorithm~\cite{bellman1958routing,ford1956network}.
Dijkstra only visits one vertex at a time and thus is work-efficient, but it is inherently sequential.
Bellman-Ford visits all vertices in a step so it requires redundant work, but can be easily parallelized.
Many parallel SSSP algorithms integrate the idea in both algorithms, and visit a subset of unsettled vertices close to the source node.
Hence, they require less work than Bellman-Ford, and have better parallelism than Dijkstra.
These algorithms are referred to as stepping algorithms (e.g., $\Delta$-Stepping\xspace and Radius-Stepping\xspace) since they process a batch of vertices in a step. This is captured by {\textsc{LaB-PQ}}\xspace{} in the stepping algorithm framework.

We present this stepping algorithm framework in \cref{alg:framework}.
This framework requires two user-defined functions, \mf{ExtDist}\xspace and \mf{FinishCheck}\xspace.
Many SSSP algorithms can be instantiated by plugging in different \mf{ExtDist}\xspace and \mf{FinishCheck}\xspace functions (see  \cref{tab:framework}).
\cref{alg:framework} starts with associate the distance array $\delta$ to a {\textsc{LaB-PQ}}\xspace{} $Q$.
It then runs in \emp{steps}.
In each step, we process vertices with distances within a threshold~$\theta$, which is computed by \mf{ExtDist}\xspace and used as the parameter of \mf{Extract}\xspace{}.
The extracted vertices will relax their neighbors using \mf{WriteMin}\xspace (\cref{line:writemin}). If successful, we call \mf{Update}\xspace{} on the corresponding neighbor. Some algorithms (e.g., $\Delta$-Stepping\xspace) contain substeps in each step.
This is captured by \mf{FinishCheck}\xspace{}---if the condition is not true, the threshold $\theta$ will not be recomputed.
We say a vertex $v$ is \defn{settled} the last time it is extracted from the {\textsc{LaB-PQ}}\xspace{} and relaxes all its neighbors (and thus its distance does not change thereafter).
We define the \defn{frontier} as all vertices in $Q$, which are those waiting to be explored to relax their neighbors.

The stepping algorithm framework applies to various algorithms as shown in \cref{tab:framework}.
We now briefly introduce them.

\myparagraph{Dijkstra and Bellman-Ford}.
Dijkstra's algorithm visits and settles the vertex with the closest distance in the frontier.
By setting $\theta$ as $\min_{v\in Q}(\delta[v])$, \cref{alg:framework} works the same as Dijkstra's algorithm, with the exception that multiple vertices with the same distances will be processed together, which does not affect correctness and efficiency.
Finding the closest vertex can be supported using $\mf{Reduce}\xspace()$ and taking $\min$ on keys.
Bellman-Ford visits all vertices in the frontier in each step, so we set $\theta$ as infinity, and in each step \cref{alg:framework} relaxes the neighbors of all vertices in~$Q$.

\myparagraph{$\Delta$-Stepping\xspace}.
As a hybrid of Dijkstra and Bellman-Ford,
$\Delta$-Stepping\xspace{} visits and settles all the vertices with shortest-path distances between $i\Delta$ and $(i+1)\Delta$ in step $i$.
Within each step, the algorithm runs Bellman-Ford as substeps.
Hence we can set $\theta$ to $i\Delta$, and use \mf{FinishCheck}\xspace to check if any newly relaxed vertex still has distance within $i\Delta$.
If not, we increment $i$ and proceed to the next step.

\myparagraph{$\Delta^*$-Stepping\xspace}. We note that \mf{FinishCheck}\xspace is not necessary for $\Delta$-Stepping\xspace, just like other stepping algorithms.
In fact, all existing implementations~\cite{zhang2020optimizing,nguyen2013lightweight,beamer2015gap,dhulipala2017} relaxed \mf{FinishCheck}\xspace in different ways. In this paper, we show that removing \mf{FinishCheck}\xspace in $\Delta$-Stepping\xspace (referred to as \emp{$\Delta^*$-Stepping\xspace}) can lead to better bounds (\cref{thm:deltass}) and good practical performance (\cref{sec:exp}).

\myparagraph{Radius-Stepping\xspace}.
In Radius-Stepping\xspace, we precompute $r_\rho(v)$, the distance from each vertex $v$ to the $\rho$-th closest vertex, for all vertices.
Then in each step, Radius-Stepping\xspace sets the threshold $\theta$ as $\min_{v\in Q}(\delta[v]+r_\rho(v))$, and then uses Bellman-Ford as substeps to compute the distances for vertices no more than the threshold.
\mf{FinishCheck}\xspace is needed by the theoretical analysis, which bounds the number of total substeps to be $O({(k_\rho n/\rho)}\cdot \log{\rho L})$.

To implement Radius-Stepping\xspace in our framework, we need an augmented {\textsc{LaB-PQ}}\xspace.
We set $r_\rho(u)$ of a vertex $u$ as the value of each record\xspace{}.
We map each record to $k+v$ for a record\xspace{} with key $k$ (distance) and value $v$ (vertex radius), and set the operator $\oplus$ as $\min$.
The threshold in \mf{Extract}\xspace is $\theta=\min_{v\in Q}(\delta[v]+r_\rho(v))$, computed by $Q.\mf{Reduce}\xspace()$.
In \cref{sec:bdpq}, we show that maintaining the augmented values does not affect the asymptotical cost bounds.

\myparagraph{$\rho$-Stepping\xspace}.
In this paper, we propose a new algorithm $\rho$-Stepping\xspace{} in the stepping algorithm framework.
$\rho$-Stepping\xspace{} extracts the $\rho$ nearest vertices in the frontier, and relaxes their neighbors.
The threshold $\theta$ is the $\rho$-th smallest element in $Q$.
We overload the notation of $\rho$ from Radius-Stepping\xspace{} because they share high-level similarities in the theoretical analysis.
The only step for $\rho$-Stepping\xspace in addition to the stepping algorithm framework is finding the $\rho$-th closest distance among all vertices in the frontier (the \mf{ExtDist}\xspace).
In our implementation, we simply use a sampling scheme that randomly pick $s=O(n/\rho+\log n)$ elements, sort them and pick the $(\rho s/n)$-th one.
More details on how to find the $\rho$-th element is in
\ifx0\undefined
the full paper,
\else
\cref{app:cpam},
\fi
and an efficient implementation is in \cref{sec:impl-theta-sampling}.

Picking the a subset of vertices with closest distances and relaxing their neighbors is not a groundbreaking idea, and has been used in the literature (e.g.,~\cite{bingmann2015bulk,zhou2019practical,alistarh2015spraylist}).
However, the extracting process in previous work is either sequential or concurrent, so none of the existing algorithms support non-trivial work and span bounds, or practical efficiency as compared to $\Delta$-Stepping\xspace.
In this paper, we argue that this simple solution can achieve both theoretical and practical efficiency.
Theoretically, we show that:
\begin{theorem}[Cost for $\rho$-Stepping\xspace]
\label{thm:rho-steps}
  On a $(k_{\rho},\rho)$-graph $G$, the $\rho$-Stepping\xspace{} algorithm has in $O\left({k_{n}m}\log\frac{n^2}{m\rho}\right)$ work and $O\left(\frac{k_n n\log n}{\rho}\right)$ span.
  If $G$ is undirected, the span is $O\left(\frac{k_\rho n\log n}{\rho}\right)$.
\end{theorem}

We will first show implementations of {\textsc{LaB-PQ}}\xspace and the cost, and then formally prove this result in \cref{sec:discussion}.
$\rho$-Stepping\xspace{} also has good practical performance, which is shown in \cref{sec:exp}.

\section{Implementation Details}\label{sec:implementation}
\label{sec:impl-bdpq}

We implemented three algorithms in the stepping algorithm framework: $\rho$-Stepping\xspace{}, $\Delta^*$-Stepping\xspace{}, and Bellman-Ford, all using array-based {\textsc{LaB-PQ}}\xspace{}.
Our implementations are simple, and are unified for the three algorithms (we only need to change \mf{ExtDist}\xspace and \mf{FinishCheck}\xspace accordingly, as shown in \cref{tab:framework}).
We present some useful optimizations we used in our implementation. Most of them apply to all the three algorithms.
Our code is available at: \url{https://github.com/ucrparlay/Parallel-SSSP}

\myparagraph{Sparse-dense optimization.}
We use sparse-dense optimization similar to Ligra~\cite{ShunB2013}.
When the current frontier is small (sparse mode), we explicitly maintain an array of vertices as the frontier.
Otherwise (dense mode), we use an array of $n$ bit flags to indicate whether each vertex is in the current frontier, and skip those not in the frontier when processing them.
The dense mode has a more cache-friendly access pattern, and avoids explicitly maintaining the frontier array, but always needs $O(n)$ time to check all vertices.
Hence, the sparse mode is used when the frontier size is small than a certain threshold.

\hide{\xiaojun{This is not used in the last version}
Our additional optimization is that we can switch from the sparse mode to the dense mode within a step. Once we detect that the next frontier size exceeds a fixed threshold, we start to only mark the bit flags instead of scattering the vertices to the next frontier (see details in the next paragraph). This saves much work in the last sparse step before dense steps.
}

\myparagraph{Queue size estimation and scattering.}
One challenge in the sparse mode is maintaining the frontier array since the size can change dramatically during the execution.
Some existing implementations (e.g., Ligra) use a parallel pack to generate the next frontier sequence, which scans all edges incident the current frontier for two rounds (one for computing offsets and another round to pack). This can incur a large overhead.
To avoid this, we use a resizable hash table to maintain the next frontier, and scatter the vertices to the next frontier by putting them into random slots in the hash table. In the process of our algorithm, we use sampling to estimate the next frontier size in order to resize the hash table.

\myparagraph{Bidirectional relaxation for undirected graphs.}
We use a novel optimization for undirected graphs.
Before the algorithm relaxes all $v$'s neighbors (\cref{line:innerloop2} in \cref{alg:framework}), it first attempts to relax $v$ using all its neighbors.
This aims to update $v$'s distance first, so it will be more ``effective'' when $v$ relaxes other vertices later.
Another reason is that parallel SSSP implementation is usually I/O bounded.
Since in relaxations, we need to check $v$'s neighbors' distances anyway, we can load them to the cache and use them to relax $v$'s tentative distance first with small cost.
This optimization only applies to undirected graphs.

\myparagraph{Threshold estimation for $\rho$-Stepping\xspace.}
\label{sec:impl-theta-sampling}
In both $\Delta$-Stepping\xspace{} and Radius-Stepping\xspace{} (although we did not implement Radius-Stepping\xspace{}), the distance threshold can be directly computed.
In $\rho$-Stepping\xspace{}, we need to compute the threshold (the $\rho$-th smallest element in the frontier) in each step.
We use the sampling-based idea as mentioned in \cref{sec:step-algo}.
In particular, at the beginning of \mf{Extract}\xspace{}, we first sample $s=O(n/\rho+\log n)$ uniformly random samples from the current frontier.
Then we sort the samples and pick the threshold from the samples.
Since $s$ is small, this step is sequential and fast.
In $\rho$-Stepping\xspace, if the frontier size is smaller than $\rho$, we pick $\theta$ as the maximum distances in the frontier.

In our experiments, we observe that in $\rho$-Stepping\xspace{}, the threshold estimation in the first dense rounds is usually inaccurate. This is because in the early stage, the $\rho$-th closest distance in the frontier is usually far from the source, and during the relaxation, much more vertices go below this threshold.
Hence, we add a heuristic to adjust the threshold: using 10\% of $\rho$ at the first two dense rounds.

\myparagraph{Large neighbor sets.}
On road networks and the begin and end for all graphs, the frontier and its neighborhood are very small.
Relaxing the neighbors in a round-based manner leads to insufficient workload and thus overhead in the synchronization cost.
To optimize this case, we use a similar ``bucket fusion'' optimization proposed by Zhang et al.~\cite{zhang2020optimizing}, which is later integrated to GAPBS~\cite{beamer2015gap}.
In our $\Delta^{*}$-Stepping\xspace and $\rho$-Stepping\xspace, when processing a vertex $v$, instead of using $v$'s direct neighbors, we run a local BFS until we reach $t=4096$ vertices (or when the tentative distances reach more than $\theta$). We use these vertices as $v$'s neighborhood $\mathcal{N}(v)$, and update them all.
Note that this information is maintained and processed locally.
As such, we can extend multiple hops in one round.
We apply this optimization in sparse rounds with average edge degree fewer than 20, and thus call them \emph{super sparse} rounds.
This optimization can greatly optimize the performance for road networks, since as shown in \cref{fig:k-rho}, the values of $k_\rho$ on road graphs is large.
The impact on the performance for social and web graphs is smaller since the dense rounds spend the most time.

\section{Introduction}
\label{sec:intro}

Given a weighted graph $G=(V,E,w)$ with $n=|V|$ vertices, $m=|E|$ edges, edge weight function $w:E\to \mathbb{R}^+$, and a source $s\in V$, the single-source shortest-path (SSSP) problem is to find the shortest paths from $s$ to all other vertices in the graph.
In this paper, we consider general positive edge weights.
Sequentially, the best known bound for SSSP is $O(m+n\log n)$ using Dijkstra's algorithm~\cite{dijkstra1959} with Fibonacci heap~\cite{fredman1987fibonacci}.
However, SSSP is notoriously hard in parallel.
Despite dozens of papers and implementations over the past decades, all existing solutions have some unsatisfactory aspects.

Practically, most existing parallel SSSP implementations ~\cite{zhang2020optimizing,dhulipala2017,nguyen2013lightweight,beamer2015gap} are based on $\Delta$-Stepping\xspace~\cite{meyer2003delta},
which is a hybrid of Dijkstra's algorithm \cite{dijkstra1959} and the Bellman-Ford algorithm~\cite{bellman1958routing,ford1956network}.
It determines the correct shortest distances in increments of $\Delta$, in step~$i$ settling down all the vertices with distances in $[i\Delta, (i+1)\Delta]$.
Within each step, the algorithm runs Bellman-Ford as substeps.

Although $\Delta$-Stepping\xspace{} is the state-of-the-art practical parallel SSSP algorithm, two challenges still remain. Theoretically, $\Delta$-Stepping\xspace{} has been analyzed on random graphs~\cite{meyer2001single,crauser1998parallelization}, but no bounds has been shown for the general case.
Practically, the parameter~$\Delta$ can largely affect the algorithm's performance.
The best choice of $\Delta$ depends on the graph structure, weight distribution, and the implementation itself.
\cref{fig:intro:delta} shows the running time of three state-of-the-art $\Delta$-Stepping\xspace{} implementations \cite{zhang2018graphit,dhulipala2017,nguyen2013lightweight} and our own $\Delta^*$-Stepping\xspace{} (a variant of $\Delta$-Stepping\xspace{}, see \cref{sec:framework}) with different $\Delta$ values, on real-world graphs (more details in \cref{sec:exp}).
A badly-chosen $\Delta$ can greatly affect the performance, and the best choices of $\Delta$ are very inconsistent for different graphs (even with the same weight distribution) and implementations.
Hence, in practice, one needs exhausting searches for $\Delta$ in the parameter space as preprocessing.

Theoretically, there has been a rich literature of parallel SSSP algorithms~\cite{ullman1991high,klein1997randomized,cohen1997using,Shi1999,cohen2000polylog,spencer1997time,blelloch2016parallel} with $o(nm)$ work and $o(n)$ span (critical path length). Most of these algorithms rely on adding shortcuts to achieve the bounds.
While these algorithms are inspiring, none of them have implementations or show practical advantages over $\Delta$-Stepping\xspace{} on real-world graphs.
We believe that one potential reason is the use of shortcuts and hopsets.
To achieve $O(n^{1-\epsilon})$ span, these algorithms need to add $\Omega(n^{1+\epsilon})$ shortcuts.
More shortcut edges contribute to more work, memory usage and footprint, hiding the advantages in the span improvement.

\ifx0\undefined
There has also been prior work on parallelizing the priority queue in Dijkstra's algorithm~\cite{brodal1998parallel,deo1992parallel,alistarh2015spraylist,sundell2005fast,linden2013skiplist,shavit2000skiplist,liu2012lock,henzinger2013quantitative,zhou2019practical,calciu2014adaptive,bingmann2015bulk,sanders2019sequential,sanders1998randomized,sanders2000fast}.
However, they do not provide interesting worst-case work and span bounds, or better performance than $\Delta$-Stepping\xspace{} in practice.

We summarize existing work on parallel SSSP in \cref{sec:related}.
\else
There has also been prior work on parallelizing the priority queue in Dijkstra's algorithm.
Early work on PRAM~\cite{brodal1998parallel,deo1992parallel} parallelizes priority queue operations, but the worst-case span bound is still $\Theta(n)$.
Other previous papers consider concurrent~\cite{alistarh2015spraylist,sundell2005fast,linden2013skiplist,shavit2000skiplist,liu2012lock,henzinger2013quantitative,zhou2019practical,calciu2014adaptive}, external-memory or other settings~\cite{bingmann2015bulk,sanders2019sequential,sanders1998randomized,sanders2000fast}.
They do not provide interesting worst-case work and span bounds, or better performance than $\Delta$-Stepping\xspace{} in practice.

We summarize existing work on parallel SSSP in \cref{sec:related}.
\fi

\input{teaser-page}

\myparagraph{Our approach}.
The three previous research directions on parallel SSSP (practical implementations, theoretical bounds, parallel priority queues) are mostly studied independently.
We aim to design parallel SSSP algorithms combining the advantages---as simple as those using parallel priority queues, achieving worst-case guarantees that match the existing bounds, and as fast as (or faster than)
$\Delta$-Stepping\xspace{} in practice.
Our key algorithmic insights include three components: a \emp{stepping algorithm framework}, which abstracts general ideas in some existing parallel SSSP algorithms, an abstract data type (ADT) \emp{\textbf{La}zy-\textbf{B}atched \textbf{P}riority \textbf{Q}ueue\xspace ({\textsc{LaB-PQ}}\xspace{})} with efficient implementations, which extracts the semantics of the priority queue needed by stepping algorithms, and two new stepping algorithms \emp{$\rho$-Stepping\xspace{} and $\Delta^{*}$-Stepping\xspace}, which are efficient both in theory and practice.

Our stepping algorithm framework (\cref{alg:framework}) abstracts the common idea in some existing ``stepping'' algorithms (e.g., Radius-Stepping\xspace{} \cite{blelloch2016parallel} and $\Delta$-Stepping\xspace{}~\cite{meyer2003delta}): in each \emph{step}, the algorithm relaxes all vertices with tentative distances within
a certain threshold, as a batch and in parallel.
The two extreme cases are the two textbook algorithms: Dijkstra's algorithm with batch size 1, and Bellman-Ford algorithm with batch size $n$.
We formalize several algorithms in this framework (\cref{tab:framework}).
Interestingly, some variants of parallel Dijkstra~\cite{bingmann2015bulk,zhou2019practical,alistarh2015spraylist} also use a similar high-level idea

The proposed ADT {\textsc{LaB-PQ}}\xspace{} abstracts the priority queue needed by the stepping algorithms.
It supports \mf{Update}\xspace{} to commit an update to the data structure, which can be lazily batched and executed in parallel. It also supports \mf{Extract}\xspace{} to return all record\xspace{s} with keys within a certain threshold in parallel.
The {\textsc{LaB-PQ}}\xspace{} is inspired by the recent work on \emph{batch-dynamic data structures}~~\cite{shun2014phase,acar2019parallel,anderson2020work,tseng2019batch,Blelloch2016justjoin,sun2018pam}, where multiple updates or queries are applied to the data structure in batches in parallel.
One advantage of {\textsc{LaB-PQ}}\xspace{} is that we do not explicitly generate the batches, but do it \emph{lazily}.
On top of the ADT, all stepping algorithms can easily use {\textsc{LaB-PQ}}\xspace's interface as a black box.
Underneath it, we provide efficient data structures to support {\textsc{LaB-PQ}}\xspace{}. We show a theoretically efficient implementation of {\textsc{LaB-PQ}}\xspace{} based on the tournament tree (\cref{sec:tour-tree}). It improves the cost bounds for existing parallel SSSP algorithms such as Radius-Stepping~\cite{blelloch2016parallel} and Shi-Spencer~\cite{Shi1999}.
In practice, we show simple implementations based on flat arrays, which makes our stepping algorithms outperform state-of-the-art software~\cite{zhang2020optimizing,dhulipala2017,nguyen2013lightweight,beamer2015gap}.

Based on the stepping algorithm framework and {\textsc{LaB-PQ}}\xspace{}, we also propose a new parallel SSSP algorithm, referred to as \emph{$\rho$-Stepping\xspace}, which is simple and efficient both in theory and in practice.
The high-level idea of $\rho$-Stepping\xspace is to relax a fixed number of unsettled vertices with small tentative distances in each step.
While a similar (but not the same) idea have been used in some parallel Dijkstra's algorithms~\cite{bingmann2015bulk,zhou2019practical,alistarh2015spraylist}, none of them have interesting bounds or practical performance comparable to $\Delta$-Stepping\xspace.
In this paper, we formally analyze $\rho$-Stepping\xspace and show work and span bounds.
$\rho$-Stepping\xspace{} achieves a better span bound than Radius-Stepping\xspace with a slightly higher work bound (\cref{thm:rho-steps}).
The work bound also applies to directed graphs (the bounds for Radius-Stepping\xspace only holds for undirected graphs).
Practically, our $\rho$-Stepping\xspace is 1.3-2.6$\times$ faster than previous implementations on social and web graphs, and is competitive on road graphs (\cref{fig:exp:heat}).

In addition to theoretical guarantees and practical performance,
another advantage of $\rho$-Stepping\xspace{} is that, it needs no preprocessing (e.g., adding shortcuts in Radius-Stepping\xspace{}) or time-consuming parameter searching (e.g., finding best $\Delta$ in $\Delta$-Stepping\xspace{}).
Our experiments (\cref{fig:rho-time}) show that, the best choice of $\rho$ is consistent and insensitive across the real-world graphs we tested

Inspired by the stepping algorithms and {\textsc{LaB-PQ}}\xspace{}, we also show $\Delta^{*}$-Stepping\xspace{}, a variant of $\Delta$-Stepping\xspace, which is simple, has non-trivial worst-case bounds (\cref{tab:bounds}), and fast in practice (\cref{fig:exp:heat}).

\hide{
\yan{The following paragraph needs revision}
$\rho$-Stepping\xspace achieves comparable or better performance as $\Delta$-Stepping\xspace{}, and similar work-span tradeoff in theory compared to Radius-Stepping\xspace{}, while also avoids expensive preprocessing or parameter searching in $\Delta$-Stepping\xspace or Radius-Stepping\xspace.
Compared to Radius-Stepping\xspace, $\rho$-Stepping\xspace has the same span bound and is only off by a little in the work bound, and is more practical.
Compared to $\Delta$-Stepping\xspace, $\rho$-Stepping\xspace is preprocessing-free, faster in almost all tested graphs, and has better theoretical guarantees on undirected graphs.
}
\myparagraph{Our Contributions.}
Combining our {\textsc{LaB-PQ}}\xspace{} with existing algorithms and our new algorithms, we achieve new bounds and efficient implementations for parallel SSSP.
These results are due to the abstraction of stepping algorithm framework and {\textsc{LaB-PQ}}\xspace{}, which greatly simplifies algorithm design, analysis, and implementation.

In theory, we show new bounds for Radius-Stepping\xspace{}~\cite{blelloch2016parallel}, Shi-Spencer~\cite{Shi1999}, $\Delta^*$-Stepping\xspace{}, and $\rho$-Stepping\xspace{}.
We note that, with no shortcuts or hopsets, it seems unlikely to show $o(n)$ worst-case span (consider a chain).
However, tighter bounds can depend on certain graph parameters, which may exhibit a good property on real-world graphs.
For example, although parallel Bellman-Ford has worst-case span of $\tilde{O}(n)$, a more precise bound is $\tilde{O}(d)$, where $d$ is the shortest-path tree depth.
Indeed, on social networks with small $d$, parallel Bellman-Ford is reasonably fast (Table \ref{tab:alltime}).
To capture this, Blelloch et al.~\cite{blelloch2016parallel} proposed a graph invariant $(k,\rho)$-graph\xspace{} that indicates how ``parallel'' a graph is.
Intuitively, a graph is a $(k,\rho)$-graph\xspace{} if every vertex reaches $\rho$ nearest vertices in $k$ hops.
We extend this concept to analyze multiple stepping algorithms.
Our experiments show that the real-world social or web graphs we tested are $(\log n,O(\sqrt{n}))$-graphs, and road graphs we tested are $(\sqrt{n},O(n))$-graphs (\cref{fig:k-rho}).
Under our framework, the stepping algorithms share common subroutines in analyses, such as the extraction lemma (\cref{lem:num-ext}) and the distribution lemma (\cref{lem:distribute}).

In practice, our framework and array-based {\textsc{LaB-PQ}}\xspace{} give unified implementations for Bellman-Ford, $\Delta^*$-Stepping\xspace{} and $\rho$-Stepping\xspace{}.
Our implementations achieve the best performance on all graphs (see \cref{fig:exp:heat}).
On the social and web graphs, $\rho$-Stepping\xspace{} is 1.3-2.6$\times$ faster than existing implementations.
On road graphs, our $\Delta^*$-Stepping\xspace is consistently the fastest and $\rho$-Stepping\xspace is competitive to previous ones.
This indicates the effectiveness of our framework since all optimizations are easily applicable to all algorithms.
We also provide an in-depth experimental study based on our framework, especially to understand the tradeoff between work and parallelism. We show how different stepping algorithms explore the frontier in steps (\cref{fig:visted-per-step,fig:enqueue}), the parameter space (\cref{fig:intro:delta,fig:rho-time}), and eventually draw interesting conclusions in \cref{sec:exp-sum}.

\hide{

The reason that we can achieve simplicity and both theoretical and practical efficiency for $\rho$-Stepping\xspace is due to the algorithmic framework and the abstraction of {\textsc{LaB-PQ}}\xspace{}.
Unlike the parallel priority queue in previous work, {\textsc{LaB-PQ}}\xspace is an ADT, not a specific data structure.
Theoretically, we can totally ignore how {\textsc{LaB-PQ}}\xspace is implemented, which significantly simplifies the $\rho$-Stepping\xspace algorithm.
Hence, we can focus on the similarities of the algorithms in the stepping algorithm framework, borrow the analysis in the recent Radius-Stepping\xspace paper~\cite{blelloch2016parallel} and apply to $\rho$-Stepping\xspace.
The practical efficiency of $\rho$-Stepping\xspace relies on the implementation of {\textsc{LaB-PQ}}\xspace, and we can borrow the highly-tuned implementation for $\Delta$-Stepping\xspace~\cite{dhulipala2017,zhang2020optimizing} and apply to our implementation.
As a proof-of-concept, the {\textsc{LaB-PQ}}\xspace data type isolates and bridges the theoretical work and practical optimizations for parallel SSSP.
}

\ifx0\undefined
Due to page limit, we postpone some analysis and full experimental results to the full version of the paper \cite{ssspfull}.
\else
\fi
We summarize our contributions of this paper as follows.

\begin{itemize}[topsep=1.5pt, partopsep=0pt,leftmargin=*
\setlength{\itemsep}{0pt}
  \item A stepping algorithm framework, which unifies multiple parallel SSSP algorithms.
  \item A new ADT {\textsc{LaB-PQ}}\xspace{} and two implementations, which are used in our analysis and implementations, respectively.
  \item A new parallel SSSP algorithm $\rho$-Stepping\xspace{}, which is preprocessing-free, simple and efficient both in theory and in practice.
  \item A new variant of $\Delta$-Stepping\xspace ($\Delta^{*}$-Stepping\xspace), which is simple, with theoretical guarantee, and fast in practice.

  \item New analyses for stepping algorithms based on $(k,\rho)$-graph\xspace{}, which include parameterized work and span bounds for $\rho$-Stepping\xspace{} (\cref{thm:rho-steps}) and $\Delta^*$-Stepping\xspace{} (\cref{thm:deltass}), and improved work bounds for Radius-Stepping\xspace{} (\cref{cor:radius}) and Shi-Spencer (\cref{cor:shispencer}).
  \item Efficient parallel implementations of Bellman-Ford, $\Delta^*$-Stepping\xspace{} and $\rho$-Stepping\xspace{}, which outperform existing ones (\cref{tab:alltime})
  \item In-depth experimental study of parallel SSSP algorithms
\end{itemize}

\hide{
The contribution is two-fold.
The first contribution is the new abstraction of \emph{batch-dynamic priority queues (BDPQ)} that bridges the different priority queue implementations and supports a variety of parallel SSSP algorithms.
This abstraction and algorithm framework are introduced in Section~\ref{sec:framework}.
Using a theoretically efficient algorithm for BDPQ introduced in Section~\ref{sec:winning-tree}, we can get improved bounds for existing algorithms.
Meanwhile, a practically efficient implementation for BDPQ gives unified and efficient implementations for $\Delta$-Stepping\xspace and Bellman-Ford.

Based on BDPQ, the second and main contribution of this paper is a new parallel SSSP algorithm \emph{$\rho$-Stepping\xspace}.
$\rho$-Stepping\xspace has several advantages: it is simple, preprocessing-free, practically efficient, and has theoretical guarantees.
$\rho$-Stepping\xspace is very simple on top of BDPQ (which can be treated as a black box).
$\rho$-Stepping\xspace is similar to practically efficient $\Delta$-Stepping\xspace and theoretically efficient Radius-Stepping\xspace, but $\rho$-Stepping\xspace does not require the expensive preprocessing step in $\Delta$-Stepping\xspace and Radius-Stepping\xspace.
Compared to Radius-Stepping\xspace, $\rho$-Stepping\xspace has the same span bound and is only off by a little in the work bound, but instead, $\rho$-Stepping\xspace is practical.
Compared to $\Delta$-Stepping\xspace, $\rho$-Stepping\xspace is preprocessing-free, almost faster in all graph instances, and has better theoretical guarantees on undirected graphs.}

\section*{Acknowledgement}
This work is in partial supported by NSF grant CCF-2103483.
\balance
\bibliographystyle{abbrv}
\input{sssp.bbl}

\ifx0\undefined
\else

\section{Preliminaries}\label{sec:prelim}

\myparagraph{Computational Model}.
{{\ifx0\undefined}}
We use the
work-span{} model for fork-join parallelism with binary forking to analyze parallel
algorithms~\cite{CLRS,blelloch2020optimal}, which is recently used in many papers on parallel algorithms~\cite{agrawal2014batching,blelloch2010low,BCGRCK08,BG04,Blelloch1998,blelloch1999pipelining,BlellochFiGi11,BST12,Cole17,CRSB13,BBFGGMS16,dinh2016extending,chowdhury2017provably,blelloch2018geometry,dhulipala2020semi,BBFGGMS18,Dhulipala2018,blelloch2020randomized,gu2021parallel}.
\else
This paper uses the work-span{} model for fork-join parallelism with binary forking to analyze parallel
algorithms~\cite{CLRS,blelloch2020optimal}, which is used in many recent papers on parallel algorithms~\cite{agrawal2014batching,Acar02,blelloch2010low,BCGRCK08,BG04,Blelloch1998,blelloch1999pipelining,BlellochFiGi11,BST12,Cole17,CRSB13,BBFGGMS16,dinh2016extending,chowdhury2017provably,blelloch2018geometry,dhulipala2020semi,BBFGGMS18,Dhulipala2018,blelloch2020randomized,gu2021parallel}.
\input{fig-overview}
\fi
We assume a set of thread{}s that share a common memory.  Each thread{} supports standard RAM instructions, and a \texttt{fork}{}
instruction that forks two new child thread{}s.
When a thread{} performs a \texttt{fork}{}, the two child thread{}s all start by running the next instruction, and the original thread{} is suspended until all children terminate.
A computation starts with a single {root} thread{} and finishes when that root thread{} finishes.
An algorithm's \defn{work} is the total number of instructions and
the \defn{span} (depth) is the length of the longest sequence of dependent instructions in the computation.
We can execute the computation using a randomized work-stealing scheduler in practice.
We assume unit-cost atomic operation \mf{WriteMin}\xspace{}$(p,\mathit{v})$\footnote{a more practical assumption is to charge $O(t)$ work and $O(\log t)$ span when $t$ operations priority update to a memory location.  It does not change the overall bound since forking $t$ parallel tasks requires $\Omega(\log t)$ span, which is already captured.}, which reads the memory location pointed to by $p$, and write value $v$ to it if $v$ is smaller than the current value.
We also use atomic operation \TAS{}$(p)$, which reads and attempts to set the boolean value pointed to by $p$ to \emph{true}{}.
It returns \emph{true}{} if successful and \emph{false}{} otherwise.

\myparagraph{Graph Notations}.
We consider a weighted graph $G=(V, E, w)$. WLOG, we assume $G$ is a connected, simple graph, with minimum edge weight $\min_{e\in E} w(e)=1$, and
no parallel edges)
We use $L=\max_{e\in E} w(e)$.
For $v \in V$, define $N(v)=\{u\,|\,(v,u)\in E\}$ as the \defn{neighbor set} of $v$.
We use $d(u,v)$ as the shortest-path distance in $G$ between two vertices $u$ and
$v$.
A \defn{shortest-path tree} rooted at vertex $u$ is a spanning tree $T$ of $G$ such that the path distance in $T$ from $u$ to any other $v\in V$ is $d(u,v)$.

\myparagraph{$(k,\rho)$-graph\xspace}. We use the concept of $(k,\rho)$-graph\xspace{} in \cite{blelloch2016parallel} to analyze stepping algorithms.
$(k,\rho)$-graph\xspace is a graph invariant highly related to the analysis of parallel SSSP algorithms.
Intuitively, a graph is a $(k,\rho)$-graph\xspace{} if any vertex can reach its $\rho$ nearest neighbors in $k$ hops.
More formally, we define the \emp{hop distance} $\hat{d}(u,v)$ from a vertex $v$ to $u$ as the number of edges on the shortest (weighted) path from $v$ to $u$ using the fewest edges.
Let $r_{\rho}(v)$ be the $\rho$-th closest distance from $v$, and $\bar{r}_k(v)$ the shortest distance from $v$ to another vertex more than $k$-hops away.

\begin{compactdef}[~$(k,\rho)$-graph\xspace~\cite{blelloch2016parallel}~]
We say a graph $G=(V,E,w)$ is a $(k,\rho)$-graph\xspace{} if for all $v\in V$, $r_{\rho}(v)\le \bar{r}_k(v)$.
\end{compactdef}

For a given graph $G=(V,E)$, we denote $k^{G}_\rho$ to be the smallest value for $k$ to make $G$ a $(k,\rho)$-graph\xspace. With clear context, we omit the superscription. $k_n$ is the shortest-path tree depth.

\myparagraph{Others.}
We use $\log n$ as a short form of $1+\log_2(n+1)$. We say $O(f(n))$ \defn{with high probability} (\whp{}) to indicate $O(cf(n))$ with probability at least $1-n^{-c}$ for $c \geq 1$, where $n$ is the input size.

\section{Related Work on Parallel SSSP}\label{sec:related}

\myparagraph{Practical parallel SSSP implementations}.
There have been dozens of practical implementations of parallel SSSP.
In this paper, we compared to a few of them.
Galois~\cite{nguyen2013lightweight} uses an approximate priority queue ordered by integer metric with NUMA-optimization
to improve the performance of SSSP.
GraphIt~\cite{zhang2018graphit, zhang2020optimizing} proposed a priority queue abstraction and a new optimization, bucket fusion, to reduce the synchronization overhead of $\Delta^*$-Stepping\xspace{}.
The optimizations are later adopted by GAPBS~\cite{beamer2015gap}, which is the one we compared to.
Julienne~\cite{dhulipala2017} proposed and used the bucketing data structure to order the vertices for
$\Delta$-Stepping\xspace{} based on semisorting~\cite{gu2015top}.
Ligra~\cite{ShunB2013} includes one of the most efficient Bellman-Ford implementations.

There has also been a significant amount of work on other implementations, include those on the distributed setting~\cite{malewicz2010pregel,zhu2016gemini,besta2017push}, GPUs~\cite{davidson2014work,wang2016gunrock}, among many others.
Our reported running time in this paper is much faster than in these papers on the same graphs, and comparing the superiorities of different settings on parallel SSSP is out of the scope of this paper.
Parallel SSSP based on parallel priority queues are reviewed in \cref{sec:bdpq}.

\myparagraph{Theoretical work on parallel SSSP}.
There has been a rich literature of theoretical parallel SSSP algorithms.
Among them, many algorithms~\cite{ullman1991high,klein1997randomized,cohen1997using,Shi1999,cohen2000polylog,spencer1997time} achieve very similar bounds to Radius-Stepping\xspace~\cite{blelloch2016parallel} we discussed in this paper, but require adding shortcut edges.
Basically the product of work and span is $\tilde{\Theta}(nm)$ (referred to as the transitive closure bottleneck~\cite{KarpR90}).
Some algorithms are analyzed based on edge weights~\cite{meyer2002buckets,meyer2001heaps}, and many others are on approximate shortest-paths~\cite{elkin2019hopsets,miller2015improved,cao2020efficient,li2020faster,andoni2020parallel} and other models~\cite{augustine2020shortest,kuhn2020computing,ghaffari2018improved}.
While these algorithms are insightful, to the best of our knowledge, none of them have implementations.

\section{Analysis for Stepping Algorithms}\label{sec:sssp}
\input{table-bounds}

With the stepping algorithm framework (\cref{alg:framework}) and {\textsc{LaB-PQ}}\xspace's implementation, we can now formally analyze the cost bounds for the stepping algorithms, which are summarized in \cref{tab:bounds}.
Our new bounds are parameterized by the definition of $(k,\rho)$-graph\xspace shown in~\cite{blelloch2016parallel}.
We first show some useful results for all stepping algorithms in \cref{sec:cost-prelim}, and use them to prove the results in \cref{sec:cost-algo}.
We later show the span for $\rho$-Stepping\xspace{} on undirected graphs in \cref{sec:undirected}, and compare with existing algorithms in \cref{sec:discussion}.

\subsection{Useful Results for All Stepping Algorithms}\label{sec:cost-prelim}

We first show two useful lemmas for all stepping algorithms.

\begin{lemma}[Number of extractions]\label{lem:num-ext}
In a stepping algorithm, a vertex $v\in V$ will not be extracted from the priority queue (\cref{line:extract} in \cref{alg:framework}) more than $k_n$ times.
\end{lemma}
\begin{proof}
  Consider the shortest path $P=\{v_0=s, v_1, v_2, \ldots, v_l=v\}$ from the source $s$ to $v$ with fewest hops.
  Since we assume the edge weights are positive, we know that $d(s,v_i)<d(s,v_{j})$ for $i<j$.
  Hence, whenever $v$ is extracted from the priority queue, the earliest unsettled vertex $v_i$ in $P$ must also be extracted and settled.
  This is because $v_{i-1}$ is already settled and have relaxed $v_i$ in previous rounds, and $d(s,v_i)\le d(s,v)$.
  Based on the definition of the $(k,\rho)$-graph\xspace, we have $l\le k_n$, which proves the lemma.
\end{proof}

\begin{lemma}[Distribution]\label{lem:distribute}
If a stepping algorithm has $S$ steps, and incurs $U$ updates (relaxations), the total work is $\displaystyle O\left(U\log(nS/U)\right)$ using tournament-tree-based {\textsc{LaB-PQ}}\xspace.
\end{lemma}

\begin{proof}
  The work of a stepping algorithm consists of modification work for relaxations (updates) and extraction work applied to the {\textsc{LaB-PQ}}\xspace{}.
  Each extracted vertex corresponds to a previous successful relaxation, and an update and an extraction have the same cost per vertex.
  Hence, we only need to analyze modification costs since extraction costs are asymptotically bounded.

  The $U$ updates are distributed in $S$ steps. Let $u_i$ be the number of relaxations applied in step $i$ ($\sum_i{u_i}=U$).
  The overall work across all steps is $W=O(\sum_i{u_i\log(n/u_i)})$.

  Since $u_i\log(n/u_i)$ is concave, $\sum_i{u_i\log(n/u_i)}\le S\cdot((\sum_i{u_i}/S)\log (n/(\sum_i{u_i}/S)))=U\log(nS/U)$, which proves the lemma.
\end{proof}

\subsection{Cost Bounds for Stepping Algorithms}\label{sec:cost-algo}

With \cref{lem:distribute,lem:num-ext} for the stepping algorithms and \cref{thm:main-tour-tree,thm:array-based} for {\textsc{LaB-PQ}}\xspace's cost, we can now show the cost bounds shown in \cref{tab:bounds} except for one given in \cref{sec:undirected}.

Dijkstra's algorithm has $O(n)$ steps and $O(m)$ relaxations, \cref{lem:distribute} gives $O(m\log(n^2/m))$ work which is essentially better than Brodal et al.'s algorithm~\cite{brodal1998parallel} (their span is also $O(n\log n)$ on the fork-join model).
Bellman-Ford has $O(k_n)$ steps and $O(k_nm)$ relaxations, so the work is $O(k_nm)$ and the span is $O(k_n\log n)$.
The following theorem shows the number of steps for $\rho$-Stepping\xspace.

\begin{theorem}[Number of steps for $\rho$-Stepping\xspace]
\label{thm:rho-steps-general}
  On a $(k_{\rho},\rho)$-graph, the $\rho$-Stepping\xspace{} algorithm finishes in $O\left(k_{n}n/\rho\right)$ steps.
\end{theorem}
\begin{proof}
  In $\rho$-Stepping\xspace{}, each step can either be a \emph{full-extract}, where $|Q|\ge \rho$ so we extract $\rho$ vertices with closest tentative distances, or a \emph{partial-extract}, where $|Q|<\rho$ so we extract all but fewer than $\rho$ vertices.
  There can be at most $O\left(k_{n}n/\rho\right)$ full-extracts, since \cref{lem:num-ext} shows that each vertex can only be extracted for $k_n$ times.
  Given that we have $n$ vertices in total, there can be at most $O\left(k_{n}n/\rho\right)$ full-extracts.
  We now show that at most $k_n$ partial-extracts can occur.
  Similar to the analysis for \cref{lem:num-ext}, once a partial-extract occurs, at least one vertex on the shortest path $P$ from source $s$ to any vertex $v$ is settled.
  Based on the definition of the $(k,\rho)$-graph\xspace, we have $|P|\le k_n$, so in total, at most $k_n$ partial-extracts can occur.
  Putting both parts together proves the theorem.
\end{proof}
Combining the result with \cref{lem:distribute} gives the work bound of $\rho$-Stepping\xspace{} in \cref{tab:bounds}.

We now show that we can get better work bounds for Radius-Stepping\xspace{} using {\textsc{LaB-PQ}}\xspace.
Radius-Stepping\xspace extracts all vertices with distance within $\min_{v\in Q}(\delta[v]+r_\rho(v))$ in each step. The original papers uses a search tree to support this operation.
We note that our {\textsc{LaB-PQ}}\xspace fully captures the need in Radius-Stepping\xspace{}.
By replacing the search tree with our tournament tree\xspace and plugging in the numbers of relaxations and steps, we get the following results.

\begin{corollary}\label{cor:radius}
Radius-Stepping\xspace~\cite{blelloch2016parallel} uses $O\left(k_\rho m\log\frac{n^2\log{\rho L}}{m\rho}\right)$ work and $O\left(\frac{k_\rho n}{\rho}\cdot \log{\rho L}\log n\right)$ span, with $O(m\log n+n\rho^2)$ work and $O(\rho\log\rho+\log n)$ span for preprocessing.
\end{corollary}

We can also improve another parallel SSSP algorithm Shi-Spencer~\cite{Shi1999} by replacing their original search-tree-based priority queue with our tournament tree\xspace{} (more details in
\ifx0\undefined
the full paper).
\else
\cref{app:shi-spencer}).
\fi

\begin{corollary}\label{cor:shispencer}
Shi-Spencer algorithm~\cite{Shi1999} can be computed using $O\left((m+n\rho)\log\frac{n^2}{m+n\rho}\right)$ work and $O\left(\frac{n\log n}{\rho}\right)$ span, with $O(m+n\rho^2\log n\log\rho)$ work and $O(\log n\log\rho)$ span for preprocessing.
\end{corollary}

We also derive the bounds for $\rho$-Stepping\xspace on directed graphs in \cref{thm:rho-steps}, and give the formal analysis for $\Delta^{*}$-Stepping\xspace{}:
\begin{theorem}\label{thm:deltass}
  $\Delta^*$-Stepping\xspace uses $O\left(\frac{k_n(\Delta+L)}{\Delta}\right)$ steps, and thus has $O\left(k_nm\log\frac{nL}{m\Delta}\right)$ work and $O\left(\frac{k_n (\Delta+L)}{\Delta}\log n\right)$ span based on {\textsc{LaB-PQ}}\xspace{}.
\end{theorem}
\ifx0\undefined
Due to the space limit, the proof of $\Delta^{*}$-Stepping\xspace is given in the full paper.
We note that for the original $\Delta$-Stepping\xspace, such bounds do not hold.
An additional factor of $k_n$ will be introduced in span if each step needs to settle down all vertices in the distance threshold.
\else
\begin{proof}[Proof of \cref{thm:deltass}]
  We first show the number of steps in $\Delta^*$-Stepping\xspace.
  The farthest vertices in the shortest-path tree has distance no more than $k_n L$ where $L$ is the largest edge weight.
  After $\lceil k_n L/\Delta\rceil$ steps, all vertices are included in the threshold, so the algorithm will finish in no more than another $k_n$ steps (the number of steps of Bellman-Ford).
  Hence, in total, $\Delta^*$-Stepping\xspace uses $S=O\left(\frac{k_n(\Delta+L)}{\Delta}\right)$ steps.

  Based on \cref{lem:num-ext}, we know the number of total relaxations are upper bounded $U=k_nm$.
  We can now use \cref{lem:distribute} to get the work bound:
  $$\displaystyle O\left(U\log\frac{nS}{U}\right)=O\left(k_nm\log\frac{nk_n(\Delta+L)}{k_nm\Delta}\right)=O\left(k_nm\log\frac{n(\Delta+L)}{m\Delta}\right)$$
  Note that since $m\ge n$, $\frac{n(\Delta+L)}{m\Delta}$ makes a difference only when $L>\Delta$.
  Hence, the work bound can be simplified as $O\left(k_nm\log\frac{nL}{m\Delta}\right)$.
\end{proof}

We note that for the original $\Delta$-Stepping\xspace, such bounds do not hold.
An additional factor of $k_n$ will be introduced if we need to settle down all vertices in each step, as shown in \cref{fig:delta-stepping-span-example}.
\fi

\subsection{Number of Steps for Undirected Graphs}\label{sec:undirected}

We can show tighter span bounds for $\rho$-Stepping\xspace on undirected graphs, which is inspired by existing results including Radius-Stepping\xspace~\cite{blelloch2016parallel} and Shi-Spencer's algorithm~\cite{Shi1999}.
\ifx0\undefined
\else
We first give the main theorem for this section.
\fi

\begin{theorem}[Number of steps, Undirected]
\label{thm:rho-steps-un}
  On an undirected $(k_{\rho},\rho)$-graph, the $\rho$-Stepping\xspace{} algorithm finishes in $O\left(k_{\rho}n/\rho\right)$ steps.
\end{theorem}

In real-world graphs, we usually have $k_\rho\ll \rho$ for large $\rho$.
Hence, by picking a large $\rho$, say $n/\log n$, $\rho$-Stepping\xspace only requires a small number of rounds and provides ample parallelism.
As a comparison, Radius-Stepping\xspace requires $O\left(\frac{k_{\rho}n}{\rho}\log\rho L\right)$ steps, a factor of $O(\log\rho L)$ more for the worst-case guarantee.

\ifx0\undefined
Due to the space limit, we show the proof in the full version of this paper, and only provide our proof sketch here.
We will show that after step $(2k_\rho+3)t$ for $t\ge 1$, $\rho$-Stepping\xspace{} will successfully settle at least the closest $t\rho$ vertices from $s$ and relax their neighbors.
This means we need $n/\rho\cdot O(k_\rho)=O\left(k_{\rho}n/\rho\right)$ steps.
We will show this by induction.
We note that the base case trivially holds when $t=1$, since $s$ can reach $\rho$ closest vertices in $k_{\rho}$ hops.
Assume this is true for $t$, we will show that this is also true for $t+1$.

More specifically, let $\mathcal{N}_{\rho}(u)$ be the set of $\rho$-nearest vertices from vertex $u$.
For simplicity, let $\mathcal{T}_\rho=\mathcal{N}_{\rho}(s)$ where $s$ is the source vertex.
The inductive hypothesis assumes that vertices in $\mathcal{T}_{t\rho}$ are settled.
We then show that within the next $2k_\rho+2$ steps, all vertices in $\mathcal{T}_{(t+1)\rho}\setminus\mathcal{T}_{t\rho}$ are settled (updated to the exact distance).

Let $v\in \mathcal{T}_{(t+1)\rho}\setminus\mathcal{T}_{t\rho}$ and $P=\{s=v_0, v_1, v_2, \dots, v_l=v\}$ be a the shortest path from $s$ to $v$ with the fewest hops.
Assume all vertices from $v_0$ to $v_i$ are in $\mathcal{T}_{t\rho}$, and beyond $v_i$ all vertices are not in $\mathcal{T}_{t\rho}$.
We show that $v$ is within $2k+2$ hops from $v_i$.
More details will be given in the full version of this paper.

\else
Our proof sketch is as follows.
We will show that after step $(2k_\rho+3)t$ for $t\ge 1$, $\rho$-Stepping\xspace{} will successfully settle at least the closest $t\rho$ vertices from $s$ and relax their neighbors.
We will show this by induction.
We note that the base case trivially holds when $t=1$, since $s$ can reach $\rho$ closest vertices in $k$ hops.
Assume this is true for $t$, we will show that this is also true for $t+1$.

We start with some notations.
Let $\mathcal{N}_{\rho}(u)$ be the set of $\rho$-nearest vertices from vertex $u$.
For simplicity, let $\mathcal{T}_\rho=\mathcal{N}_{\rho}(s)$ where $s$ is the source vertex.
The inductive hypothesis assumes that vertices in $\mathcal{T}_{t\rho}$ are settled.
We now show that within the next $2k_\rho+2$ steps, all vertices in $\mathcal{T}_{(t+1)\rho}\setminus\mathcal{T}_{t\rho}$ are settled (updated to the exact distance).

Let $v\in \mathcal{T}_{(t+1)\rho}\setminus\mathcal{T}_{t\rho}$ and $P=\{s=v_0, v_1, v_2, \dots, v_l=v\}$ be a the shortest path from $s$ to $v$ with the fewest hops.
Assume all vertices $v_0$ through $v_i$ are in $\mathcal{T}_{t\rho}$, and beyond $v_i$ all vertices are not in $\mathcal{T}_{t\rho}$.
We will first show that $v$ is within $2k+2$ hops from $v_i$.

Throughout the analysis, we consider $v_j$ on path $P$ where $j=i+k+2$ as a special vertex, if it exists. If not, it directly means that $v$ is no more than $k+2$ hops away from $v_i$.
We use the neighbor set of $v_j$ to show $l\le i+2k_\rho+3$ and complete the proof.
To start, we show the following lemma.
  \begin{lemma}
  \label{lem:nonn1}
  $v_{i+1}$ is not in $\mathcal{N}_{\rho}(v_j)$.
  \end{lemma}
  \begin{proof}
  From the definition of $(k,\rho)$-graph\xspace, all vertices in $\mathcal{N}_{\rho}(v_j)$ are within $k$ hops from $v_j$.
  Since the shortest path with fewest hops from $v_{i+1}$ to $v_j$ has $k+1$ hops, $v_{i+1}$ cannot be in $\mathcal{N}_{\rho}(v_j)$.
  \end{proof}

We now show the following lemma that $\mathcal{T}_{t\rho} \cap \mathcal{N}_{\rho}(v_j)=\varnothing$.
It says that no vertices in $\mathcal{T}_{t\rho}$ are close enough to be processed in the previous steps.
We use \cref{lem:nonn1} and $v_{i+1}$ as a separating vertex.

  \begin{lemma}
  \label{lem:nonn2}
  None of the vertices in $\mathcal{T}_{t\rho}$ is in $\mathcal{N}_{\rho}(v_j)$.
  \end{lemma}
  \begin{proof}
    Assume to the contrary that there exists a vertex $u\in \mathcal{T}_{t\rho}\cap \mathcal{N}_{\rho}(v_j)$.
    Then we know $d(s,u)<d(s,v_{i+1})$, since $u$ is one of the $t\rho$ closest vertices from $s$ ($u\in \mathcal{T}_{t\rho}$) but $v_{i+1}$ is not.
    From \cref{lem:nonn1}, we know $v_{i+1}\notin \mathcal{N}_{\rho}(v_j)$.
    This means that if $u\in \mathcal{N}_{\rho}(v_j)$, then $d(u,v_j)<d(v_{i+1}, v_j)$.
    Combining the above two conclusions, we know that $d(s,u)+d(u,v_j)<d(s,v_{i+1})+d(v_{i+1},v_j)$, which leads to a contradiction that $v_{i+1}$ is on the shortest path from $s$ to $v_j$.
  \end{proof}

  \cref{lem:nonn2} says none of the $\rho$ closest vertices of $v_j$ are in $\mathcal{T}_{t\rho}$.
  We next show that for all vertices in $\mathcal{T}_{(t+1)\rho}\setminus\mathcal{T}_{t\rho}$ are close to $v_j$.

  \begin{lemma}
  For any $u\in \mathcal{T}_{(t+1)\rho}\setminus\mathcal{T}_{t\rho}$, it is either $v_{i+1}$, or it is within $k$ hops from $v_j$.
  \end{lemma}
  \begin{proof}
    Again, assume to the contrary that $u$ is not within $k$ hops from $v_j$ and is not $v_{i+1}$.
    In that case, $u$ is not in $\mathcal{N}_{\rho}(v_j)$.
    In that case, any vertex $u'\in \mathcal{N}_{\rho}(v_j)$ should be closer to $s$ than $u$.
    Since $u\in \mathcal{T}_{(t+1)\rho}$, any $u'\in \mathcal{N}_{\rho}(v_j)$ should also be in $\mathcal{T}_{(t+1)\rho}$.
    Lemma \ref{lem:nonn2} shows that none of $\rho$ vertices in $\mathcal{N}_{\rho}(v_j)$ is in $\mathcal{T}_{t\rho}$. Therefore, all the $\rho$ vertices in $\mathcal{N}_{\rho}(v_j)$ and $u$ should be in $\mathcal{T}_{(t+1)\rho}\setminus\mathcal{T}_{t\rho}$, which indicates that $|\mathcal{T}_{(t+1)\rho}\setminus\mathcal{T}_{t\rho}|>\rho$, leading to a contradiction.
  \end{proof}

As a result, we know that $v\in \mathcal{T}_{(t+1)\rho}\setminus\mathcal{T}_{t\rho}$ is at most $2k+2$ hops away from a vertex in $\mathcal{T}_{t\rho}$, i.e., $l\le i+2k+2$.
Recall that $j=i+k+2$, so within $2k+2$ hops from $v_i$, we can get the shortest distance of any $v\in \mathcal{T}_{(t+1)\rho}\setminus\mathcal{T}_{t\rho}$.

  \begin{lemma}
  \label{lem:updaterho}
  After step $(2k_\rho+3)t+h+1$, vertex $v_{i+h}$ must have been settled and have relaxed all its neighbors for $h\le (l-i)$.
  \end{lemma}
  \begin{proof}
  The inductive hypothesis indicates that $v_i$ must have relaxed all its neighbors in the first $(2k_\rho+3)t$ steps. Therefore, as its neighbor, $v_{i+1}$ should be settled at step $(2k_\rho+3)t+1$.

  Next, we show that $v_{i+1}$ will be extracted from the {\textsc{LaB-PQ}}\xspace{} to update all its neighbors no later than step $(2k_\rho+3)t+2$.
  We know that $v_{i+1}$ is no farther than $v$ from $s$, and $v$ is in $\mathcal{T}_{(t+1)\rho}\setminus\mathcal{T}_{t\rho}$.
  This means that there cannot be at least $\rho$ unsettled vertices in the frontier that are closer to $s$ than $v_{i+1}$, so $v_{i+1}$ will be extracted.
  Similarly, $v_{i+2}$ will be settled in step $(2k_\rho+3)t+2$. We can show this inductively that the lemma holds for all $h\le l-i$.
  \end{proof}
  Plugging in $h=l-i$, we can know that $v=v_l$ must have been settled and relaxed all its neighbors before step $(t+1)(2k+3)$.
  With \cref{lem:updaterho}, we directly get \cref{thm:rho-steps-un}.

Lastly, we note that \cref{thm:rho-steps-un} is an upper bound.
\cref{lem:updaterho} shows that $v_{i+h}$ is settled no later than step $(2k_\rho+3)t+h+1$, but it can be settled earlier.
This is shown in \cref{fig:visted-per-step} by our experiment, and on real-world graphs, the number of steps is very small.
\fi