Parallel and (Nearly) Work-Efficient Dynamic Programming
Xiangyun Ding Yan Gu Yihan Sun
















































  University of California, Riverside
  

xding047@ucr.edu



  University of California, Riverside
  

ygu@cs.ucr.edu



  University of California, Riverside
  

yihans@cs.ucr.edu


























The idea of dynamic programming (DP), proposed by
Bellman in the 1950s, is one of the most important algorithmic
techniques. 


However, in parallel, many fundamental and sequentially simple problems become
more challenging, and open to a (nearly) work-efficient solution (i.e.,
the work is off by at most a polylogarithmic factor over the best sequential solution). 
In fact, sequential DP algorithms employ many advanced optimizations such as
decision monotonicity or special data structures, and achieve
better work than straightforward solutions. 
Many such optimizations are inherently sequential, 
which creates extra challenges for a parallel algorithm 
to achieve the same work bound. 

The goal of this paper
is to achieve (nearly) work-efficient parallel
DP algorithms by parallelizing classic, highly-optimized and
practical sequential algorithms.
We show a general framework called the  for parallel DP algorithms,
and use it to solve several classic problems. 
Our selection of problems includes Longest Increasing Subsequence (LIS), sparse Longest Common Subsequence (LCS), 
convex/concave generalized Least Weight Subsequence (LWS), Optimal Alphabetic Tree (OAT), and more. 

We show how the  can be used to achieve the same level of optimization as the sequential algorithms,
and achieve good parallelism. 
Many of our algorithms are conceptually simple, and we show some experimental results as proofs-of-concept. 































[500]Theory of computation Parallel algorithms



















 
 


















Parallel and (Nearly) Work-Efficient Dynamic Programming


    Yihan Sun
    Received March XX, 2024; accepted XX XX, 2024
============================================================







§ INTRODUCTION


The idea of dynamic programming (DP), since proposed by Richard Bellman in the 1950s <cit.>, has been extensively used in algorithm design, 



and is one of the most important algorithmic techniques. 
It is covered in classic textbooks and basic algorithm classes, and is widely used in research and industry. 
The goal of this paper is to achieve (nearly) work-efficient (defined below) and parallel DP algorithms 
based on parallelizing classic, highly-optimized and practical sequential algorithms.

At a high level, a DP algorithm computes the DP values for a set of states (labeled by integers) by a recurrence.



The recurrence specifies a set of transitions from state j to state i, i.e., how D[j] can be used to compute D[i][


More generally, a transition may compute D[i] from multiple other states. 
All algorithms in this paper only requires one state j in the transition to compute D[i].
].

We call j a decision at i. 
D[i] is then computed by taking the best (minimum or maximum) among all decisions,
which we call the best decision at i.

Throughout the paper, we will use i^* to denote the best decision of state i. 
We introduce more concepts about DP in <ref>. 


At a high level, a DP algorithm solves an optimization
problem by breaking it down to simpler subproblems, memoizing the answers to the subproblems, and using them to find the answer to the original problem. 
The subproblems are usually indexed by integers, referred to states. 
The DP value of a state is determined either by a boundary condition (i.e., initial values),
or from other states, specified by a DP recurrence. 
The recurrence specifies a set of transitions from state j to state i, i.e., how D[j] can be used to compute D[i].
We call j a decision at i. 
D[i] is then computed by examining all decisions at state i, and taking the best (minimum or maximum) among them. 
We call the decision that gives the final value of D[i] the best decision at i, and denote it as i^*. 


One can view the states and transitions as a directed acyclic graph (DAG), which we refer to as a DP DAG. 
In this DAG, each vertex is a state, and an edge j to i denotes a transition from j to i. 
Since such an edge indicates that computing D[i] logically requires D[j], 
we also call it a dependency, and say i depends on j. 
Sequentially, we can compute all states based on a topological ordering. 
For simplicity, we always assume that the (integer) order of the states is a valid topological ordering.


Unfortunately, many DP algorithms (even some simple ones sequentially) are hard to parallelize,
and are especially hard to achieve work-efficiency (the work asymptotically matches the best sequential algorithm) or even near work-efficiency 
(off by a polylogarithmic factor). 
We note that on today's multicore machines with tens to hundreds of processors, 
achieving low work is one of the most crucial objectives for designing practical parallel algorithms. 

One particularly intriguing and somewhat ironic challenge for achieving work-efficient parallel DP algorithms is that, 
sequential algorithms are extremely well-optimized. 
In many cases, an optimized DP algorithm does not need to process all edges (transitions) in the DP DAG; 
some even do not need to process all vertices (states).
We review the literature at the end of this section. 
In fact, almost all textbook DP solutions can be optimized to achieve lower work than the straightforward solution 
that directly computes the DP values of all states based on the recurrence. 
Such examples include longest increasing subsequence (LIS), (sparse) longest common subsequence (LCS), 
(convex/concave) least weight subsequence (LWS), and many others discussed in this paper. 












First, sequential DP algorithms usually process all states iteratively. 
When processing a state sequentially, it is easy to guarantee that all states it depends on have been finalized. 
Therefore, each state can be computed exactly once to obtain its true DP value. 
Despite the iterative essence of the sequential DP algorithm, 
it is easy to observe that a state does not need to wait for all states before it to finish, but only those it depends on. 
This allows for much shallower dependency structure depth than trivially O(n), where n is the number of states. 
Therefore, many existing ideas achieve parallelism by
identifying as many ready (i.e., all their decisions are finalized) states as possible and process them together. 
Since this essentially breaks the sequential order of the states,
an important challenge for work-efficiency is to avoid processing a state before it becomes ready, or to
carefully bound such “wasted” work on touching unready states. 








As mentioned, such optimizations of interest usually allow the algorithm to skip certain edges and vertices in the DP DAG.
In the literature, some optimizations that can avoid processing all states are called sparsification <cit.>.
In this paper, we call it vertex sparsification since it skips vertices in the DP DAG. 
One typical example is the LCS problem. 
In a nutshell, although the DP recurrence contains O(mn) states (m and n are the input string lengths), 

only a subset of them need to be computed. 
We will discuss LCS as an motivation example in <ref>. 

The second type of optimization is to strategically skip some edges, i.e., avoid using certain transitions. 
Following the literature on (vertex) sparsification, we call these optimizations edge sparsification. 
One typical example is the optimized LIS algorithm (e.g., the one in TAOCP <cit.>)—
when processing a state i, 
the algorithm maintains a data structure to precisely find the best decision of i in O(log n) cost.
Effectively, only O(n) transitions are processed in the algorithm, requiring O(nlog n) work. 




One typical DP optimization that is both theoretically elegant and practically useful 
is decision monotonicity (DM).
At a high level, DM indicates that,
if a state i has best decision i^*,
then another state j>i must have its best decision j^*≥ i^*, called convex case, 
or the concave case[The definitions of convexity and concavity are interchanged in some other papers.]
where j^*≤ i^* or j^*≥ i.
At a high level, DM indicates that two states i and j>i must have their best decisions j^*≥ i^*, called the convex case[The definitions of convexity and concavity are interchanged in some other papers.],
or the concave case where either j^*≤ i^* or j^*≥ i (see <ref>). 




Hence, when finding the best decision for state j,
one can narrow down the possible range of j^* based on the known best decisions of previous states, 
and thus avoid processing all transitions. 
DM has been widely studied in the sequential setting (e.g., <cit.>),
and is also closely related to concepts such as 
quadrangle inequalities <cit.> and Monge property <cit.>.
Sequentially, 
using DM saves a polynomial factor than the naïve DP algorithm in many applications <cit.>. 
In the parallel setting, however, among the papers we know of <cit.>,
most from the 90s, very few of them take advantage of DM to reduce work. 
Indeed, none of them are work-efficient, and most of them have a polynomial overhead, 
which limits their potential applicability on today's multicore machines. 
The only nearly work-efficient results <cit.> focus on the concave case of one specific problem.


In parallel, however, 
among the papers we know of concerning DP with DM optimizations (e.g., on edit distance <cit.>, xxx), 

none of them are (even nearly) work-efficient when comparing to the sequential solutions using DM. 


The challenge of using DM in parallel lies in two aspects. 
First, sequentially we skip the transitions for state i by observing the best decisions of all states before i.
When multiple states are processed in parallel, they cannot see each other's best decision, 
making it hard to skip the same set of “useless” transitions as in the sequential setting.
Second, many classic sequential DP algorithms with DM relies on efficient data structures such as monotonic queues, which are inherently sequential. 
Achieving the same work bound in parallel also requires careful redesign of the underlying data structures. 




In this paper, we study parallel DP algorithms to achieve the same work as highly-optimized sequential algorithms. 

Given a sequential algorithm Γ with certain optimizations, 
our goal for work-efficiency is to process (asymptotically) the same number of transitions and states as in Γ. 
Regarding parallelism, we hope to achieve the best possible parallelism indicated by the transitions processed by Γ.

We formalize our goals in <ref>. 


Our solution is based on an algorithmic framework that generally applies to almost all DP algorithms.
We call this framework the . 

At a high level, our framework specifies how to correctly identify a subset of states
that do not depend on each other and process them in parallel. 

We then present how to do so efficiently for specific problems.
To achieve work-efficiency, our key ideas are two-fold. 
First, many of our algorithms use prefix-doubling 
to bound the additional work on processing unnecessary states.
Second, we design new parallel data structures to skip unnecessary transitions. 




This paper studies general approaches for parallel DP, 
with a special focus on applying the non-trivial, effective optimizations found in the sequential context to parallel algorithms. 
We select classic DP problems and their optimized sequential solutions, and parallelize them using novel techniques. 
Our framework unifies one existing parallel LIS algorithm <cit.>, 
and provides new parallel algorithms for various problems such as sparse LCS,
convex/concave generalized LWS and GAP edit distance (GAP), and optimal alphabetic tree (OAT).
All the algorithms are (nearly) work-efficient with non-trivial parallelism. 
Among them, we highlight our contributions on parallelizing the DP algorithms with decision monotonicity. 
The core of our idea is a parallel algorithm for convex/concave generalized LWS. 
We apply it to other problems such as GAP and OAT, and achieve new theoretical results. 
For OAT, we partially solve the open problem in <cit.> by providing 
a work-efficient algorithm with polylogarithmic span with input as positive integers with word size n^(n). 
We present our theoretical results in <ref>. 
We believe this is the first paper that achieves near work-efficiency in parallel for a class of DP algorithms with DM. 


































Although the main focus of this paper is to achieve low work in theory, 
an additional goal is to make the algorithms simple and practical. 
We implement two algorithms as proofs-of-concept (code available at <cit.>). 
On 10^9 input size, 
both of them outperform sequential solutions 
when the depth of the DP DAG is within 10^5,
and achieve 20–30× speedup with smaller depth of the DP DAG. 



Due to page limit, some proofs are provided in the full version of the paper <cit.> (see the supplementary material). 




Related Work.
Dynamic programming (DP) is one of the most studied topics in algorithm design.
The seminal survey by <cit.> reviewed two types of optimization techniques sequentially, including decision monotonicity (e.g., <cit.>) and sparsity (e.g., <cit.>).
In this paper, we mainly focus on parallelizing the sequential algorithms in this scope, and we will review the literature of each problem in the corresponding section.

DP is also widely studied in parallel.

There exists rich literature on optimizing various goals in different models, 
such as span (time) in PRAM (e.g. <cit.>),
I/O cost in the external-memory/ideal-cache model (e.g., <cit.>), rounds in the MPC model <cit.>, or on the BSP model <cit.>. 
However, these papers either only considered the DP algorithms without the optimizations, 
or incur polynomial work overhead, except for <cit.> for one specific problem.
Instead, our paper tries to parallelize the efficient and practical sequential algorithms while maintaining low work. 
Some other works try to parallelize certain types of DP or applications using DP (e.g., <cit.>).
Alternately, our work aims to provide a general approach to parallelize almost all DP algorithms.

Notations. 
We use O(f(n)) with high probability () in terms of n as O(cf(n)) with a probability of at least 1-n^-c for c≥ 1. 

We use Õ(·) to hide polylog(n) where n is the input size. 




Given a list of states D_i and a weight function w(i,j), in most cases, a dynamic programming algorithm computes the value of D_j by taking the minimum (or maximum) of D_i+w(i,j) from the feasible region of state D_j. 
Such a relationship among the states is called a DP recurrence.
We refer to i as the decision of the state D_j if D_j=D_i+w(i,j).




A series of seminal work on dynamic programming is decision monotonicity (DM)[Properties with similar ideas include quadrangle inequalities and Monge property.] (cite them all), and we will review the literature in more details in the next section.
The high level idea is that, the decisions for computing certain states can be used to limit the feasible regions of other states without affecting the correctness.
As such, algorithms taking advantages of DM can saves the total work to compute the recurrences, usually by Ω̃(n) for input size n.
Most of the “natural” weight functions in real-world applications such as biology, geology, and AI lead to decision monotonicity <cit.>.

The parallelism of dynamic programming has been extensively studied in the past decades.
Again in most cases, a DP computation can be considered as matrix multiplications on the (+,min)-semiring and computed in polylogarithmic span (longest dependence chain), but requires at least cubic work (number of operations) and thus is not practical.
There have been lots of studies on optimizing the span of a DP computation (cite Galil's and other's papers) and the I/O efficiency (cite Rezaul's and other's papers).
However, most of these algorithms do not consider decision monotonicity (DM).
Hence, for applications with the DM properties, these algorithms also incur higher work bounds than the best sequential counterparts that consider DM (i.e., these parallel algorithms are not work-efficient).



We are aware of a few papers on DP that considers DM, mainly focusing on edit distance (cite AALM and the other ones) and some other specific problems (cite Lawrence's papers).
Unfortunately, these algorithms are also not work-efficient when comparing to the best known sequential algorithms using DM properties on these topics.

Given the prevalence of multicore machines, work-efficiency is generally adopted as the golden rule for practical parallel algorithms.
Particularly in this paper, we want to study and understand the parallelism in the large group of (sequential) DP algorithms using decision monotonicity, without sacrificing additional work.

Many work-efficient (or nearly work-efficient) parallel DP algorithms have been recently proposed, on topics such as edit distance <cit.>, LIS and others <cit.>.
One main challenge in parallelizing DP algorithms is that the longest dependence chain in the computation, which is longest chain for the decisions, is usually linear to the input size in the worst case, leading to limited parallelism.
As such, achieving sublinear span bound sometimes needs to bring in additional work.
An alternative option is to directly parallelizes the sequential algorithms <cit.>. The goal is to achieve parallelism in the original computation, and a careful design of the parallel schemes will ideally not cause much additional work.
Such an parallel algorithm should process all objects in the proper order based on the dependences—it should 1) process as many objects as possible in parallel (as long as they do not depend on each other), and 2) process an object only when it is ready (all objects it depends on are finished), to avoid redundant work.
To formalize the two requirements, we say an algorithm is round-efficient <cit.> if its span is Õ(D) for a computation with longest logical dependence length D.
We say an algorithm is work-efficient if its work is asymptotically the same as the best sequential algorithm.




In a high level, a DP algorithm solves an optimization
problem by breaking it down to simpler subproblems, and memoizing the answers to the subproblems to find the answer to the original problem. 
The subproblems are usually labeled by integers, called the states of the original problem.
The DP value of a state is determined either by a boundary condition (i.e., initial values),
or from other states, specified by a DP recurrence. 
The recurrence computes the DP value of state i, noted as D[i],
from several options, each corresponding a previous state j. 
The final DP value will be the best among all the options. 
We call all such states j that 




 



§ MODEL AND FRAMEWORK




We use the work-span model in the classic multithreaded model with binary-forking
<cit.>.
We assume a set of threads that share the memory.
Each thread acts like a sequential RAM plus a fork instruction that forks two threads running in parallel.
When both threads finish, the original thread continues.
A parallel-for is simulated by forks in a logarithmic number of steps.
A computation can be viewed as a DAG.
The work W of a parallel algorithm is the total number of operations,
and the span (depth) S is the longest path in the DAG.
In practice, we can execute the computation with work W and span S using a randomized work-stealing scheduler <cit.> in time W/P+O(S)
with P processors with high probability.
A parallel algorithm is work-efficient, if its work is O(W), where W is the work of the best known or the corresponding sequential algorithm,
and nearly work-efficient if its work is Õ(W). We use Õ(·) to hide polylog(n) where n is the input size. 




In the following, we start with presenting some basic concepts and optimizations of DP.
We then present our algorithmic framework, the , for parallel DP algorithms. 
We note that  by itself does not guarantee efficiency. 
Later in <ref> we present how to use more techniques to enable (near) work-efficiency for the
problems of interest. 




 §.§ Basic Concepts in Dynamic Programming



A dynamic programming (DP) algorithm solves an optimization problem by breaking it down to simpler subproblems, 
and memoizing the answers to the subproblems to find the answer to the original problem. 
The subproblems are usually labeled by integers, called the states of the original problem.
The algorithm computes a DP array D[·], 
where D[i] memoizes the answer to the state (subproblem) i. 
In this paper, we use i or i to denote a state i. 
We call D[i] the DP value for state i. 
The DP value D[i] is computed either by a boundary condition, 
which is an initial value, or computed by DP values of other states, specified by a DP recurrence. 
In this paper, we study recurrences in the following form. 

    D[i]=min/max_j∈(i) f_i,j(D[j])




In the recurrence, D[i] can be computed from several other states (noted as a set (i)), 
called the decisions at i. 
D[i] chooses the best among all decisions. 
If j is a decision at i, we say i depend on j. 


We say state j^* is the best decision of state i, if D[i]=f_i,j(D[j^*]), i.e.,
state j gives the optimal solution for state i. We use array [i]=j^* to record the best decision for each state i. 



A dynamic programming (DP) algorithm, as introduced in <ref>, 
usually involves computing the DP value D[i] of each state i based on the recurrence, starting from the boundary conditions. 
We use i to denote a state i, or directly use i with clear context. 


A DP algorithm solves an optimization
problem by breaking it down to subproblems, memoizing the answers to the subproblems, and using them to find the answer to the original problem. 
The subproblems are usually indexed by integers, referred to as states. 
With clear context, we directly use i to refer to “state i”.
The DP value of a state is determined either by a boundary condition (i.e., initial values),
or from other states, specified by a DP recurrence. 
This paper studies recurrences in the following form:

    D[i]=min/max f_i,j(D[j])





where j is a decision at i.
Function f_i,j(·) indicates how the DP value of state j can be used to update state i.






The transitions (i.e., dependencies) among states form a DAG G=(V,E) as introduced in <ref>. 


We use ji to denote an edge from j to i in the DAG. 






In the recurrence, D[i] can be computed from several other states (noted as a set (i)), 
called the decisions at i, 
and chooses the best (min or max) one among them.
In our paper, we consider the scenario where each decision at i only depend on one other state. 
If i's DP value can be computed from j, we say i depends on j. 






During a DP algorithm, we may maintain the DP value of a state and update it by the recurrence. 
We call the process of updating D[i] by D[j] a relaxation, or say j relaxes i. 
A relaxation is successful if the DP value is updated to a better value, i.e., a lower (higher) value if the objective is minimum (maximum).


We call the actual DP value of a state the true or finalized DP value to distinguish from the DP value being updated during the algorithm, 
which we call the tentative DP value. 

We say a state i is finalized if we can ensure that its true DP value has been computed, and tentative otherwise.

Among all tentative states, we say a state is ready, if all its decisions are finalized, and unready otherwise. 

A  DP algorithm will process all transitions and states based on a topological ordering.




DP Optimizations. 
Instead of computing the recurrence straightforwardly, 


many algorithms can optimize the computation by skipping vertices and/or edges in the DAG to save work. 
We call such algorithms optimized DP algorithms. 






For example, given an input sequence A[1..n], the optimized LIS algorithm <cit.>
maintains a data structure to precisely find the best decision of each state in O(log n) cost, 
and only processes O(n) transitions instead of O(n^2) as suggested by the recurrence. 
Similarly, given two input sequences A[1..n] and B[1..m], 
the optimized LCS algorithm only needs to process all states D[i,j] where A[i]=B[j],
instead of O(nm) states in the recurrence. 
Another typical optimization is decision monotonicity where the best decisions of previous states can narrow down the range for best decisions for later states, 
which skips transitions and saves work.
Typical examples include concave/convex generalized LWS (see <ref>), OAT (see <ref>), GAP (see <ref>), etc. 

We will discuss all these algorithms in this paper. 



    
  * Precise Relaxation. Some DP algorithms can directly find the best decision of a given state once all its predecessors are finalized. 
  This is usually accomplished by certain data structures, such as a binary search structure in LIS <cit.>. 
    
  * Decision monotonicity (DM). 
        DM indicates that the best decisions of previous states can narrow down the range for potential best decisions for later states, 
        and therefore skips transitions and saves work.
        Typical examples include concave/convex generalized LWS (see <ref>), GAP (see <ref>), etc. 
  
  
  




 §.§ Commonly-Used DP Optimizations

Instead of computing the recurrence straightforwardly, 
many algorithms can optimize the computation by modifying the DAG based on some rules. 

We call such algorithms optimized algorithms. 

Such optimizations can possibly skip vertices and/or edges in the DAG to save work. 



We discuss two well-adopted DP optimizations as examples below.





  
  * Edge Refinement. Some optimizations can narrow down the scope of potential best decisions for a given state.
  
  This usually skips unnecessary edges in the DAG, i.e., 
  we can avoid evaluating edge ji if j has been excluded from the best decision of i. 
    We discuss two types of optimization in this category.
        
    
  * Precise Relaxation. The first case is to directly find the best decision of a given state once all its predecessors are finalized. 
  This is usually accomplished by certain data structures, such as a binary search structure in LIS <cit.>. 
    
  * Decision monotonicity (DM). 
    
    
    
        
        DM indicates that the best decisions of previous states can narrow down the range for potential best decisions for later states, 
        and therefore skips transitions and saves work.
        Typical examples include concave/convex generalized LWS (see <ref>), GAP (see <ref>), etc. 
        
  
  
  * Shortcutting. Some optimizations shortcut two states that are not directly connected in the original DP DAG,
  which skips both vertices and edges. 
  
  
  In many cases, the shortcuts are only conceptual and do not need to be evaluated. 
  We also discuss two types in this category, summarized as follows.
  
    
  * Greedy Choices. In some applications, among all the decisions at a state, 
  a greedy choice can be proved optimal. 
  In this case, a series of greedy decisions can be made and we can shortcut the corresponding edges in the DP DAG.
  Effectively, this skips some states (vertices) in the DAG and thus saves work. 
  Typical examples include LCS and the Vishkin-Landau edit distance algorithm. We discuss LCS as an example in <ref>.
    
  * Edge Associativity. Some other existing divide-and-conquer (DaC) DP algorithms <cit.>, although not covered by this paper, 
    can also be broadly categorized as shortcutting algorithms.
    
    Unlike the previous case where a series of greedy decisions from j to i are directly employed, 
    these algorithms shortcut j to i by enumerating and finding the best path between them.
    The idea applies to associate state transitions, i.e., we can combine the transition from j to l and that from l to i to 
    compute D[i] from D[j]. 
    
    
    This is usually achieved by
    enumerating possible intermediate states l, typically situated midway on the path between j and i. 
    It then recursively finds the optimal solutions from j to l and l to i and combine the solutions, forming a DaC strategy. 
    This usually leads to polylogarithmic depth of the DP DAG (and thus the span bound). 
    
    
    
    
    
    
    
  
  
  






We refer to the first category as Edge Sparsification.
Usually in these applications, |(i)| for i is large (e.g., O(n) size).
These optimizations aim at reducing the size of |(i)| (usually to O(1) or O(log n)).
Commonly used techniques include:



  
  * Precise Relaxation. For some recurrences, it is possible to use some efficient data structures to directly and precisely find the best decision of a given state. A typical example is to use range queries, as in the longest increasing subsequences (LIS) and activity selection problem <cit.>add more. 
  
  * Decision Monotonicity. Different from the previous technique, decision monotonicity narrows down the range of [i] by using best decisions of the other states. It is generally more challenging to parallelize since the search range of each state is not known ahead of time and determined on-the-fly during the sequential execution.
  
  






  
  * Vertex Sparsification. In some applications, among all the decisions at a state, 
  a greedy decision can be proved to be optimal. In this case, a series of greedy decisions can be made and we can short-circuiting the corresponding edges in the DP DAG.
  Effectively, this skips some states (vertices) in the DAG and thus saves work. 
  Typical examples includes LCS and the Vishkin-Landau edit distance algorithm. We discuss LCS as an example in <ref>.
  
  
  
  * Edge Sparsification. In certain applications, it is possible to narrow down the scope of potential best decisions for a given state.
  Therefore, we can avoid evaluating edge ji if j has been excluded as the best decision of i. 
  This paper studies two types of edge sparsification, summarized as follows.
  
    
  * The first case is to directly and precisely find the best decision of a given state once all its predecessors are finalized. 
  This is usually accomplished by use certain data structures, such as a binary search structure or various search trees in LIS and activity selection problem <cit.>add more. 
    
  * The second and more interesting case is the decision monotonicity (DM) property. As introduced in <ref>, this property indicates that given any state i with best decision i^*, the best decision j^* for another state j>i must satisfy j^*>i^*, which is called the convex case, or j^*>i^*, called the concave case. 
        In this paper, we focus on the convex case since convex DP recurrences are more widely observed in real world applications. 
        Typical examples include 1D cluster problem (defined in <ref>), 2D GAP problem (defined in <ref>), etc. 
        DM indicates that the best decisions of previous states can help to narrow down the range for potential best decisions for later states. 
        By doing this, one can skip edges corresponding to the unpromising decisions and save work. 
  
  
  Some other existing work, although not covered by this paper, can also be broadly categorized as combinations of vertex and edge sparsification. 
  Typical examples include divide-and-conquer algorithms such as AALM edit distance algorithm and the 
  and retrieve a path by enumerating the vertex on a boundary and concatenating the two paths in both subproblems.  Typical solutions include (cite AALM, the other, and those for LIS).
  
  
  
  







Although not our main focus, we would also introduce the other category of optimizations, which we refer to as Vertex Sparsification.  
The high-level idea is to avoid computing all states in the original recurrence, which will inevitably modifies the edges to guarantee correctness.
Typical examples are when the number of states is large but |(i)| is small (e.g., longest common subsequences), and commonly seen approaches include:


  
  * Shortcutting. Some algorithms add shortcut edges which combine a series of dependencies and replace them. A typical use case is to shortcut a list of 0-weight edges, such as the Landau-Vishkin edit distance algorithm <cit.>.
  
  * Short-circuiting. Another common approach is to enable divide-and-conquer, and retrieve a path by enumerating the vertex on a boundary and concatenating the two paths in both subproblems.  Typical solutions include (cite AALM, the other, and those for LIS).



Note that these techniques do not guarantee fewer edges, but they usually reduce the depth of a DAG, and can be used to parallelize certain sequential algorithms.
Also, the two categories are not disparate.  
For example, in <ref>, LIS with O(n) vertices and O(n^2) edges can equivalently be solved by LCS with O(n^2) vertices and edges.
Hence, the classic O(nlog n) algorithm is an edge sparsification (precise relaxation) on LIS but a vertex sparsitication on LCS.
In this paper, we focus on Decision Monotonicity with is an edge sparsification technique.





This paper mainly parallelizes sequential algorithms based on decision monotonicity and greedy choices. 
We also show an example of applying our framework to precise relaxation in <ref>, which directly gives an existing algorithm <cit.>.






 §.§ Parallelizing Sequential DP Algorithms


We now discuss our goal to parallelize a sequential DP algorithm.
Our primary goal is to achieve (asymptotically) the same computation as an optimized sequential algorithm. 


More formally, for a sequential algorithm Γ computing a recurrence R with certain optimizations, 
we define the 
Γ-optimized DAG, 
denoted as G_Γ=(V_Γ,E_Γ), as follows. 
G_Γ is the same as the DP DAG on R with some edges highlighted: 
for all edges that are processed by Γ, 
we highlight them in G_Γ, and call them the effective edges. 
We call the other edges normal edges.





These effective edges can be used to fully complete the computation. 

For a sequential DP algorithm Γ, we hope its faithful and best possible parallelization Λ to:


  

  
  * process the same effective edges in G_Γ and achieve the same work as Γ, and
  
 
  
  * have span proportional to the effective depth of G_Γ, 
  defined as the largest number of effective edges in any path in G_Γ. 



Namely, if our goal is to parallelize Γ and perform the same computation, 
the parallel algorithm Λ has to process all edges ji processed by Γ. 
This means i can be finalized only after j is finalized, 

so the span is related to the effective depth as defined above. 


We use E_Γ to denote the set of effective edges,  

and (G_Γ) as the effective depth of G_Γ. 







More formally, we define the following concepts.

We say a parallel algorithm Λ is an optimal parallelization of a sequential DP algorithm Γ, if 
1) the set of edges processed in Λ and Γ satisfies E_Γ⊆E_Λ and |E_Λ|=O(|E_Γ|), 
2) the work of the two algorithms are asymptotically the same, and
3) the span of Λ is Õ((G_Γ)). 
Namely, Λ can process a superset of edges of that in Γ, but not asymptotically more,
and its span is proportional to the effective depth of the Γ-optimized DAG.


In addition, we define the perfect parallelization.
In an omniscient version of Γ, we only need to process the best decisions based on their dependencies. 
We define the Γ-perfect DAG, denoted as G^*_Γ, 
as subgraph of G_Γ that only contains the best decision edges (both effective and normal edges).




We say an optimal parallelization Λ of Γ is also a perfect parallelization 
if the span of Λ is Õ((G^*_Γ)). 






In addition, we define the perfect parallelization of a sequential algorithm Γ as follows.



In an omniscient version of Γ, we only need to evaluate the best decisions based on their dependencies. 
We define the perfect depth of a state s as follows.
Starting from s, we will chase the best decision path in the DAG until we see a state initialized by the boundary condition.
The number of effective edges on this path is the perfect depth of s.
The perfect depth of a DAG G, denoted as ^* is the largest perfect depth among all states.
An optimal parallelization Λ of Γ is also a perfect parallelization 
if the the span of Λ is Õ(^*(G_Γ)). 


We say a perfect parallelization Λ is a perfect parallelization of Γ, if 

1) Λ is an optimal parallelization of Γ, and
2) the span of Λ is Õ(^*(G_Γ)). 


While the definitions seem abstract, we will later show that they are intuitive for concrete problems. 
For example, in longest common subsequence (LCS), both (G_Γ) and (G^*_Γ) are the output LCS length k 
(but for other algorithms they can be different), 
and our goal is to achieve Õ(k) span for a perfect parallelization. 




In addition, we define the perfect parallelization of Γ as follows.


Consider an omniscient version of algorithm Γ, denoted as Γ^+, that knows the best decision of all states and only evaluates them in a proper order. 
Then the DAG G_Γ^+ should only have the edges from [i] to i as effective, 
and also have all shortcut edges (as normal edges) from j to i as long as there exists a path from j to i in the original DAG. 
We say a parallel algorithm Λ is a perfect parallelization of Γ, if 
1) the set of edges evaluated in Λ satisfies |E^*(G_Λ|)=O(|E(G^*_Γ^+)|), and
2) the span of Λ is Õ(^*(G_Γ^+)). 


In other words, the omniscient algorithm Γ^+ forsees all dependencies between states, but only evaluates the useful ones. 
Therefore, without additional knowledge such as edge association, ^*(G_Γ^+) is the best dependency depth we can achieve. 
We will show that many of our new algorithms in this paper is not only optimally parallelized, but also perfectly parallelized. 

While the definitions seem abstract, we will later show that these concepts are intuitive for concrete problems.
The effective depth ^*(G_Γ^+) suggests the longest list of best decisions that one depends on the other. 
For example, in LCS, ^*(G_Γ^+) is exactly the output LCS length k, 
and our goal is to achieve Õ(k) span for a perfect parallelization. 




In addition, we define the concept of perfect parallelization.
We consider the DAG G_Π=(V_Π, E_Π), which only consists of edges from [i] to i for a given recurrence,
with possible shortcuts. 
We say a parallel algorithm Λ is a perfect parallelization of a recurrence, if 
1) the set of edges evaluated in Λ satisfies E_Π⊆ E^*_Λ and |E(G_Λ|)=O(|E(G_Π)|), and
2) the span of Λ is Õ((G_Π)), where (G_Π) is the depth of G_Π. 

In an omniscient approach, the state i only relies on [i], and can be computed when [i] is ready.
We will show that many of our new algorithms in this paper is not only optimally parallelized, but also perfectly parallelized. 



Note that the “perfect parallelization” of a sequential algorithm does not directly suggest optimality in work or span bounds for the same problem. 
One can possibly achieve better bounds by redesigning the recurrence and/or sequential algorithm with fewer edges or a shallower depth.


Instead of finding or redesigning a different DAG to obtain new optimizations, 
our focus is to provide parallelization of existing sequential algorithms with optimizations. 







A parallel algorithm Λ is a perfect parallelization of a sequential algorithm Γ, if 
1) Λ evaluates the same set of dependencies as Γ, i.e., G_Γ=G_Λ, and
2) the span of Λ is Õ(𝑑𝑒𝑝𝑡ℎ(Γ)). 




Challenges.
We identify three major challenges in such parallelizations, particularly for DP with decision monotonicity.

The first challenge is to process multiple vertices in parallel.
Recall that in edge sparcification, we prune an edge ji if j cannot be the best decision of i.





However, when multiple states are processed in parallel,
the states in the same batch cannot see the result of each other. 
Therefore it can be hard to narrow down the range for [i] in the same way as the sequential algorithm. 

The second challenge is to check the readiness of vertices.
To achieve work efficiency, we need to accurately decide all ready vertices to proceed in each round, and process them in parallel.
Particularly for DP with decision monotonicity, the readiness can be affect by any other vertices.
How to find all vertices work-efficiently is usually highly nontrivial.

Third, many such optimizations relies on efficient data structures, such as monotonic queues. 
These data structures are inherently sequential. To enable the same work as the sequential algorithm, we have
to design new parallel data structures.

Hence, we provide a comprehensive and synergetic solution by a novel algorithmic framework (as shown in <ref>), and multiple algorithmic insights and data structures (discussed in <ref>) to tackle these challenges, and introduce them in the rest of this paper.


In the following, we first present our new algorithmic framework: the ,
which provides a correct, although not necessarily efficient parallelization for general DP algorithms. 

On top of it, for each specific problem, 
we will show how to achieve low work, 
which, as discussed in <ref>, can be highly non-trivial.




 §.§ Our Framework: the 


Our idea is based on the phase-parallel framework <cit.> (see below)
adapted to DP algorithms. 



The phase-parallel framework by Shen et al. aims to identify (as many as possible) operations that do not depend on each other, and process them in parallel. 
Directly applying this framework to DP algorithms will give the following algorithm outline:



[style=mystyle]

  
   While there exist any tentative states:
  
    
  * Find the set of ready states as 
    
  * Process all states in  in parallel and mark all of them as finalized 
    
  



We call each iteration of the while-loop a round. 
We call the set of states being processed in round i the frontier of round i, noted as _i. 
While the phase-parallel framework gives a high-level approach in achieving parallelism, it does not indicate how to do so (i.e., how to identify the ready states in each round).



We now introduce our , which uses a novel approach to identify the frontier _i in each round, particularly for a DP computation G_Γ. 


The  identifies the unready tentative states and put s on them; then it uses all s
to outline a cordon to mark the boundary of the frontier. We summarize the algorithm in the following steps.
Note that every step can be processed in parallel. 


  

  
  *  Mark all states as tentative and initialize them by the boundary condition.
  

  
  *  
  If a tentative state j can successfully relax another tentative state i (i.e., update D[i] to a better value), 
  put a  on state i. Such a  means that all the descendants (inclusive) of state i are unready. 
  We say this state and all its descendants are blocked by the  in this case. 
  
      Therefore, a state is ready if there is no  on any of its ancestors (inclusive). 
  

  
  *  For each ready state, use its DP value to relax the tentative DP values of its descendants. Usually, we need to do so implicitly to achieve efficient work. We discuss more details later.
  

  
  *  Mark all ready states as finalized. Clear all the sentinels.
  

  
  *  If there still exist tentative states, go to <ref> and repeat. 


We will first prove that the algorithm is correct, i.e., it computes correct DP values for all states. 
Later we will show some motivating examples to help understand this algorithm in <ref>. 


 



  The  is correct.


  This is equivalent to show that, when we find a ready state and later mark it finalized, its DP value must be finalized. 
  
  We will show this inductively. At the beginning, the ready states are those with zero in-degree in the DAG, and their true DP values are specified by the boundary cases. 
  Clearly, their DP values cannot be relaxed by other states and they will be identified as ready in our algorithm. 
  Their DP values are also computed by the boundary in <ref>. Therefore, our algorithm is correct in the first round. 
  
  Assume the algorithm is correct up to round r. We will show that round r+1 also correctly finds the ready states and their DP values. 
  Assume to the contrary that a state i identified in <ref> does not have its true DP value yet. 
  
  This means that the best decision of i has not relaxed i to the finalized value. 
  Let j=[i]. 
  Note that j cannot be finalized: if so, before j is marked as finalized in <ref>, in <ref> it must have relaxed i.
  We discuss two cases. 
  First, j is finalized. This is impossible because before j is marked as finalized in <ref>, in <ref> it must have relaxed i, which leads to a contradiction.
  
  
  
  Therefore, j is tentative. 
  In this case, let j_0=i, and state j_x be the best decision of state j_x-1 for x≥ 1, i.e., we chase the chain of best decisions and get a list of states i=j_0, j_1=[j_0], etc.
  Let us find the first x such that D[j_x] is not the true DP value but D[j_x+1] is the true DP value. 
  We then prove that state j_x+1 must be a tentative state. If we consider j_y as the first finalized state on this chain, 
  then j_y-1 is a tentative state, and also must already have its true DP value (because j_y has relaxed it in <ref>). 
  Note that since state j_y-1 is a tentative state, all states between j_0 and j_y-1 must be tentative. 
  Since j_y-1 has its true DP value, and i=j_0 does not have its true DP value, 
  the first switch point j_x must be between j_0 and j_y-1, and must also be a tentative state.
  
  Therefore, state j_x+1 is a tentative state that can relax state j_x, so it will put a  on state j_x. As a descendant of j_x, 
  state i must be blocked by j_x, and cannot be identified as a ready state. This leads to a contradiction.   
  Therefore, if a state i is identified as ready in our algorithm, its DP value must have been finalized.   






The  tells which states/vertices should be in the frontier in each round, but the algorithm does not show how to do so efficiently. 
Especially in <ref>,
it is almost infeasible to use the finalized DP values to explicitly update all other states.
In this case, we have to develop new parallel data structures to facilitate this step. 
Next, we first use LIS and LCS as two examples to illustrate our framework.
We then apply it to more involved cases that use decision monotonicity in <ref>.




§ MOTIVATING EXAMPLES ON LIS/LCS




To help the audience understand the more complicated algorithms using DM in the following sections,
we first provide two simple examples on , especially on how to compute the cordon efficiently.

Longest Increasing Subsequence (LIS). We first use LIS as an example.
Given an input sequence A_i, LIS computes the maximum of D[i] such that:

    D[i]=max{1, max_j<i, A_j<A_i D[j]+1}

We use n as the input size and k as the output LIS length. 
Directly computing this recurrence takes O(n^2) work.
Sequentially, we can process all states one by one, and compute max_j<i, A_j<A_i D[j] using a binary search structure, giving O(nlog k) total cost <cit.>. 
The binary search precisely finds the best decision of each state, and only n transitions are processed. 

Many parallel LIS solutions have also been proposed <cit.>.
We will show that solving LIS using our  framework will essentially
give an existing algorithm <cit.> 
and is a perfect parallelization of the sequential O(nlog k) algorithm. 

Based on , the boundary case is to set all tentative DP values as 1. 
Then we will attempt to use the current tentative DP values for relaxation. 
In this case, for a state i, as long as there is any other state j<i 
with A_j<A_i, D[i] can be relaxed to a better value 2. 
Therefore, all ready states are those input objects that are prefix-min
elements in the sequence, i.e., A_i is the smallest value among all A_1..i.
We can set these states as ready, update all other states and repeat.
Note that since all unfinalized states have the same tentative DP value of 2,
we do not need to explicitly update the values in D[i], but
can just maintain a global variable as the current tentative DP value.
By the same idea, the ready states in the next round would be the prefix-min elements
in the input after removing the finalized states. 
The same observation (repeatedly finding prefix-min elements) is exactly the core idea of the algorithm in <cit.>.
In their algorithm, they further use a tournament tree to identify prefix-min elements efficiently and achieve efficient O(nlog k) work, 
and O(klog n) span. 
Note that k is exactly the perfect depth of the DP DAG, which is the longest dependency between best decisions. 



  Combining with a tournament tree,  leads to a perfect parallelization of sequential O(nlog k) LIS algorithm in <cit.>,
  where n is the input size, and k is the LIS length.



Longest Common Subsequence (LCS).
LIS has a close relationship with other important problems such as LCS. 
Here we will revisit our cordon-based LIS algorithm from the view of LCS, which also leads to a new parallel LCS algorithm. 


Given two sequences A[1..n] and B[1..m] (m≤ n), LCS aims to find a common subsequence C of A and B such that C has the longest length among all
common subsequences. 
An LIS problem can be reduced to an LCS problem by first relabeling all input elements by 1..n based on their total order,
and then finding the LCS between this new sequence and a sequence B=⟨1,2,…, n⟩. See <ref>(a)–(c) for an illustration.
LCS has been extensively studied both sequentially <cit.>
and in parallel <cit.>.
The standard DP solution defines each state D[i,j] as the LCS for A[1..i] and B[1..j], and uses the following recurrence:

    D[i,j]=
               0,     i=0  j=0 
    
               D[i-1,j-1]+1,     A[i]=B[j] 
    max{D[i-1,j],D[i,j-1]},    .


These transitions correspond to horizontal, vertical and diagonal edges on a grid (see an example in <ref>(c)). 
A known sequential optimization (i.e., sparsification) <cit.> to this recurrence is to observe that only the edges correspond to the diagonal edges with A[i]=B[j] are useful. 
The computation is equivalent to finding the longest path from the bottom-left to top-right corresponds to these effective (red) edges,
and all other edges and states can be skipped. This can lead to a sequential algorithm with O(log n) cost,
where  is the number of pairs (i,j) such that A[i]=B[j]. For LIS, there are exactly =n such effective edges.

We will show how  can be used to parallelize this optimization. 
Starting with the boundary where D[i,j]=0, we observe that the DP value D[i,j] of any state with A[i]=B[j]
can be updated to a better value. Therefore, we will put a  at 
each of such states to indicate that they should be updated. All such s will block the top-right part of the 
grid. In this way, the blocked region is clearly marked by a staircase region, as shown in <ref>(d). 
Therefore, the entire region within the first cordon has the DP value 0.
By repeatedly doing this, we will effectively find that the region between the cordons of round i+1 and i are those
states with DP value (LCS length) i. The algorithm finishes in k rounds where k is the LCS length.

The problem boils down to efficiently identifying the cordon (i.e., the staircase) in each round.
Note that in LIS there is at most one effective edge in each column (see <ref>(d)),
while in LCS there can be multiple effective edges in each column (see <ref>(e)).
We will show an interesting modification to the original LIS algorithm that can handle this more complicated setting.
Here we sort all edges by column index as the primary key (from the smallest to largest) and row index as the secondary key (but from largest to smallest).
An example is illustrated in <ref>(f).
Then, we will still use a tournament tree to maintain this list, and apply prefix-min on the row indexes.
It is easy to see that a state/edge is on the cordon if its row index is smaller than or equal to the prefix-min.
A tournament tree can identify, mark, and remove these states in O(llog (L/l)) work and O(log n) span <cit.>, where l is the number of diagonal edges on the cordon.
We thus have the following theorem.


  Combining with a tournament tree,  leads to a perfect parallelization (O(log n) work and O(klog n) span) of sequential LCS algorithm in <cit.>,
  where n and m<n are the input sequence sizes,  is the number of pairs (i,j) such that A[i]=B[j], and k is the LCS length. 


Since L=O(n^2), the O(log L) terms in the cost of the tournament trees is stated as O(log n) in the theorem.


Interestingly, to the best of our knowledge, this is the first parallel LCS algorithm with o(mn) work and o(min(n,m)) span for sparse LCS problem (i.e., =o(mn) and k=o(min(n,m))). 
Meanwhile, this algorithm is quite simple—we provide our implementation in <cit.> and experimentally study it in <ref>.
Another interesting finding is that our algorithm implies how to map LCS to LIS (previously only the other direction is known).
Given two input strings A and B, if we sort all (i,j) pairs for A_i=B_j by increasing i (primary key) and decreasing j (secondary key), then LCS is equivalent to the LIS on the secondary keys (the j(s)) of this sorted list.

We will show more sophisticated parallelization of DP algorithms in the next sections.
Our LCS algorithm will be a subroutine in the more involved parallel GAP algorithm introduced in <ref>.






Here we will use the longest increasing subsequence (LIS) as the example, which has recently received much attention in the parallel algorithm community <cit.>.
All of these algorithms are solving it on an equivalent form on the longest common subsequence (LCS) problem as shown in <ref>(c) and discussed below.
Given an input sequence A_i, LIS computes the maximum of D[i] such that:

    D[i]=max_j<i, A_j<A_i D[j]+1

Directly computing this recurrence takes O(n^2) work, and one can accelerate it by using a range search to compute max_j<i, A_j<A_i D[j].
Sequentially, if we compute D[i] in the increasing order of i, then we only need 1D range-max queries and updates, and any balanced BST suffices here with O(nlog n) total work.
However, this solution is inherently sequential.
To parallelize it, existing solution view LIS as an equivalent LCS form, and an example shown in <ref>(a–c).
By replacing each element as a diagonal edge in a 2D-grid based on (i,A_i), computing LIS is equivalent to compute the LCS, or the longest paths from (0,0) to (n,n) regarding the diagonal edges.

The  for LIS.
The recent parallel LIS algorithm by Gu et al. <cit.> can be viewed as the  in the 2D LCS view (<ref>(c–d)).
For this application, there are n vertices of interest in this 2D grid, each is the pointed by a diagonal arrow.
Recall that in the , our goal is to decide the set of ready vertices _i such that all vertices they depend on (inferred by the <ref>) are finalized.
The geometric view here is that, all unfinished vertices will block the vertices on the top-right corner.
Hence, we put a sentinel on each unfinished vertex, and they will form a cordon that denotes the boundary between the ready and unready vertices.
For instance, in <ref>(d), the gray shaded region is the ready region for _1.
The next cordon is decided by the rest 5 vertices (sentinels), and we highlight the next cordon _2.
This cordon can be computed by taking the prefix-min for all unfinished vertices from the bottom row to the top.
Another property here is that, all vertices in _i has the DP value as i.
In total, this will trivially gives an O(nk) work and O(klog n) span algorithm where k is the LIS length.
Gu et al. <cit.> showed that by using a parallel (static) tournament tree, finding a cordon with m vertices on it can be computed in O(mlog (n/m)) work and polylogarithmic span.
The overall work of the algorithm is therefore O(nlog(n/k)) and the span is Õ(k).

Extensions to LCS.
We note that by formalizing the process of cordon finding, the above algorithmic idea can extended to a general LCS problem.
The only difference is that every diagonal in the 2D grid can exist (<ref>(e)).
However, this difference does not affect the high-level approach of cordon finding in the , but just with minor details in checking and maintaining the readiness.
Due to the space limit, we will discuss this algorithm in Appendix <ref>, and only give the theorem here.

  Given m diagonal edges in the 2D grid, LCS can be computed in O(mlog(m/k)) work and O(klog m) span, where k is the LCS length.

Interestingly, similar approaches are categorized in vertex sparsification (see <ref>) and studied sequentially in <cit.>.



Let P[i] be the smallest index of arrows in column i.
In each round, the shape of the staircase is formed by the prefix-min of P[1..n], denoted by P_min[1..n].
Then for column i, all arrows under P_min[i-1] are eliminated in this round,
and we update P[i] for the next round.
For example, in <ref>(e) we have P_min[6]=4, so in the 7'th column the staircase eliminates two arrows,
and we update P[7] to 7 as the next arrow in column 7.
In each column, if the arrows are pre-sorted, we can figure out how many arrows are eliminated by a dual binary search,
where the cost is proportional to the number of arrows eliminated.
The total number of arrows is O(k(m+n)).


  If all the arrows are sorted in each column,  leads to a perfect parallelization (O(k(m+n)) work and Õ(k) span) of sequential LCS algorithm in <cit.>, 
  where n and m are the input sequence sizes, and k is the LCS length. 

















§ PARALLEL GENERALIZED LWS


We now discuss the convex/concave generalized least weight subsequence (GLWS) problem, 
which is one of the most classic cases of decision monotonicity (DM).

Given a cost function w(j,i) for integers 0≤ j<i≤ n and D[0], 
the GLWS problem computes

    D[i]=min_0≤ j<i{E[j]+w(j,i)}

for 1≤ i ≤ n, where E[j]=f(D[j],j) can be computed in constant time. 
The original least weight subsequence (LWS) problem <cit.> is
a special case when E[j]=D[j].  
Here we use the general case E[j]=f(D[j],j) that
has the same sequential work bound <cit.>, 
because the generalized version is needed in many applications (see examples in <ref>).
The GLWS problem is also referred to as 1D/1D DP by <cit.>.



The GLWS problem is highly relevant to other important problems (e.g., line breaking <cit.>,
optimal alphabetic trees <cit.>, and a number of computational geometry problems <cit.>).

The essence of GLWS is to cluster a list of 1D objects based on spatial proximity and minimize the total weighted sum of all clusters.

Other problems such as optimal paragraph formation <cit.>
and finding optimal alphabetic trees <cit.> can be reduced to LWS problem. 
A number of computational geometry problems also reduce to 1D/1D DP <cit.>.
As an intuitive example, consider selecting a subset of villages on a road (with their coordinates known) to build post offices to minimize the total cost, 
where 
w(j,i) is the cost of using one post office to serve villages j+1 to i. 


This gives a GLWS problem with D[i] as the lowest cost to serve the first i villages and E[i]=D[i]. 
The DP recurrence enumerates all possible decisions j such that the last post office serves the villages j+1 to i,
and takes the minimum cost among all possible decisions j. 
Practical cost functions w (e.g., a fixed cost plus a linear or quadratic cost to the service range or sum of distances from villages to the post office) are convex, 
which implies DM—for two states i and j>i, their best decisions i^* and j^* satisfy j^*≥ i^*.
Symmetrically we can show that for concave cost functions w, either j^*≥ i or j^*≤ i^* holds, although concave cost functions are less common in the real-world applications of GLWS.

Many practical cost functions w, e.g., w(j,i) is the total distance from the post office to all villages in this range,
are convex (defined below) and lead to DM. 

This means that for two states i and j>i where their best decisions are i^* and j^*, then we must have j^*≥ i^*. 

One can also define a concave cost function, where we must have either j^*≥ i or j^*≤ i^*. 

Given its high relevance in practice, 

convex GLWS has been studied in parallel. 
<cit.> showed an algorithm with O(n^2log n) work and O(log^2n) span. 


Later, <cit.> showed an improved algorithm with O(n^1.5log n) work and O(√(n)log n) span.
Despite the interesting algorithmic insights in these algorithms, the polynomial overhead in work limits their potential to outperform the classic sequential solutions with Õ(n) work <cit.>.
For the concave case, some works <cit.> achieve near work-efficiency and polylog span on the original LWS,
but the ideas cannot be applied to generalized LWS.

In this section, we show how to use the  to parallelize a well-known sequential GLWS algorithm with O(nlog n) work, which works for both convex and concave DM. 
Although efficiently applying  here requires many sophisticated algorithmic techniques, 
our parallel algorithm (<ref>) remains practical and it indeed outperforms the sequential algorithm in a wide parameter range (see <ref> for details).
It is also the key building block for many other algorithms shown later in <ref>.

We start with preliminaries and the classic sequential algorithm, then discuss how to use  to parallelize it.
We will use the convex case when describing the algorithm since it is used more often in practice, and discuss the concave case in <ref>.








 §.§ Preliminaries


Convexity, Concavity and Decision Monotonicity.
The convexity of the cost function w is defined by the Monge condition <cit.>.
We say that w satisfies the convex Monge condition (also known as quadrangle inequality <cit.>) if for all a<b<c<d,

    w(a,c)+w(b,d) ≤ w(b,c)+w(a,d).

We say that w satisfies the concave Monge condition (also known as inverse quadrangle inequality) if for all a<b<c<d,

    w(a,c)+w(b,d) ≥ w(b,c)+w(a,d).




Consider two states i and j>i with best decisions i^* and j^*. 
A convex weight function leads to DM such that j^*≥ i^*. 
A concave weight function leads to DM such that either j^*≥ i or j^*≤ i^*.


Another condition closely related to the Monge condition is the total monotonicity <cit.>.
We say a n× m matrix A is convex totally monotone if for a<b and c<d,

    A(a,c)≥ A(a,d) ⇒ A(b,c)≥ A(b,d).

We say a n× m matrix A is concave totally monotone if for a<b and c<d,

    A(a,c)≤ A(a,d) ⇒ A(b,c)≤ A(b,d).

Let c_i be the column index such that A[i,c_i] is the minimum value in row i.
The convex total monotonicity of A implies that r_1≤ r_2≤...≤ r_n,
while in the concave case r_1≥ r_2≥...≥ r_n <cit.>.
Also, if A is convex totally monotone, any submatrix B of A is also convex totally monotone.
The convex/concave Monge condition is the sufficient but not necessary condition for convex/concave total monotonicity.

In GLWS, we call f_i,j=E[j]+w(j,i) a transition from j to i.
The convex/concave decision monotonicity is equivalent to the convex/concave total monotonicity of f_i,j.
Note that if w satisfies the convex/concave Monge condition, so does f_i,j.
But the convex/concave total monotonicity of w does not guarantee the convex total monotonicity of f_i,j.
Throughout this paper, we assume the convex/concave Monge condition of w,
but all our theorems only need the convex/concave total monotonicity of f_i,j.



A more general condition to enable DM is convex/concave total monotonicity <cit.>.
Our algorithm works as long as the transition function is convex/concave totally monotone. 
For page limit, we provide the full background of related discussions in <ref><cit.>.


The Sequential Algorithm.
The best (sequential) work bound for convex GLWS is O(n) <cit.>, and O(nα(n)) for the concave case <cit.>. 


However, both of them are mainly of theoretical interest since they are complicated and have large constants in both work and space usage. 




We parallelize a simpler and more practical algorithm with O(nlog n) work <cit.>.


This algorithm computes D[1..n] in order.

It implicitly maintains the best decision array [1..n].
When the algorithm finishes computing D[i], the algorithm updates [i+1..n] using D[i],
then [j] (j>i) will be the best decision of state j among states 0 to i. 





However, maintaining and updating this array of size n for n iterations require quadratic work.

Observe that after computing D[i], [i+1..n] must be non-decreasing in the convex case, 
and must be non-increasing in the concave case <cit.>. 
Hence, the algorithm maintains a “compressed” version of [i+1..n] by a list of triples ([l,r],j),
which indicates that all states between l and r have best decision j, i.e., ∀ i'∈ [l,r], [i']=j. 
The list is maintained by a monotonic queue, 

which is a classic data structure based on double-ended queue, and is inherently sequential. 
In the i-th iteration, we can directly find the decision of state i from the queue. 

After obtaining D[i], the monotonic queue can be updated in O(log n) amortized cost to consider i as the best decision for all later states.

In total, this algorithm costs O(nlog n) work.
Here we refer to the audience to the original paper for details of this algorithm. 
We will call this algorithm . 
Making use of DM,  only processes transitions between each state i and its best decision. 

The DAG G_ for this algorithm includes normal edges ji for all j<i,
and exactly n effective edges [i]i for all states i. 


To efficiently batch-update the best decision array,
the sequential algorithm maintains double-ended queue storing triples ([l,r],j),
which means that [l..r]=j.
The triples in the queue are sorted in ascending order of [l,r]
so that [i] can be read directly from the triple in the front of the queue.
After we compute D[i],
in the convex case we can start from the back of the queue (or from the front in the concave case)
and check whether i is better than the current best decision.
For a interval we only check its two endpoints, and if i can successfully relax both of them we can pop out this interval and proceed to the next one.
Otherwise, we do a binary search inside this interval to find the cutting point.



Duering the process we add O(n) intervals to the queue,
so the push/pop of the invervals taks O(1) time amortized for each i.
Plus the cost of binary search, the total work is O(nlog n).


Due to simplicity, this algorithm is usually the choice of implementation in practice. 
We will show a parallel version of this algorithm using . 


Existing Parallel Algorithms
For the LWS problem where f(D[j])=D[j] and w(j,i) is concave,
there exists a near work-efficient algorithm with polylogarithmic span <cit.>.
We are unaware of existing (nearly) work-efficient and well-parallelized algorithm for the generalized LWS, 
whether for convex or concave cases.




 §.§ Parallel Convex GLWS


 


 
 

We first give the parallel algorithm of convex GLWS. 
We will use the “post-office” problem mentioned above as a running example to explain the concepts, 
but our algorithm works for general cases. 

Following the idea of the phase-parallel algorithm, 
with the current finalized states, the goal is to find all ready states as the frontier, 
where their true DP values can be computed from the finalized ones. 
We will use our  to find the frontier in each round. 
Naïvely, the recurrence suggests that a state depends on all states before it.


However, note that a state is essentially ready as long as its best decision has been finalized. 
For the convex case, we will use the fact as shown below. 



  In convex GLWS, let S={j: j>i∧[j]≤ i}, which is the set of states with best decisions no later than 
  state i; then S is a consecutive range of states starting from i+1. 




This is a known fact in the sequential setting (can be proved by induction). It suggests that
the frontier of each round in the phase-parallel algorithm is a consecutive range of states. Based on this idea, 
we will maintain  as the last finalized state in each round. 
Then in the next round, ideally, the algorithm should find the cordon at , 
where all states [+1, -1] are ready and can compute their true DP value from (i.e., have their best decisions at) states no later than . 
We show an example of the post-office problem to illustrate the phase-parallel framework in <ref>.

Based on the discussions above, ideally, in round i, the ready states in the frontier
are those where the best solution contains i post offices. 
This is because their best decision must have i-1 post offices, and must be finalized in the previous round. 

This high-level idea is presented in the main function of <ref>. 
Starting from =0, given the current finalized states [0,], 
we will find all ready states [+1,-1] using the , 
which essentially will find the cordon at . 
We explain this part with more details in <ref>. 



Similar to the sequential algorithm,
we also maintain a data structure B to store all triples ([l,r],j) in order,
which indicates that all states between l and r have best decisions at j. 
This data structure is essential to guarantee the efficiency of finding the next frontier, 
and also has to be updated after each round with the new DP values (<ref>).






We still follow the basic idea in the sequential algorithm to maintain [1..n].
The difference is that in each round,
we need to detect all ready states so they can be computed in parallel.
<ref> shows the process.
Suppose D[1..] have been finalized.
By the , the  is the minimum index of all sentinels.
If we compute , states [+1..-1] will be ready.
After we compute D[+1..-1], we use their values to update [..m].
Note that as  depends on some state between [+1..-1],
its best decision must be updated, and so do all states after .
<ref> illustrates one round of <ref>.



Let k=^*(n) be the number of times the best[] function must be iteratively applied before the result is 0.
For <ref> we have the following conclusion:


  If w satisfies the convex Monge condition, <ref> ends in no more than k rounds.



  Consider the sequence [v_0,v_1,v_2,...,v_k] where v_0=0, v_k=n and [v_i]=v_i-1 for i∈[1,k].
   starts with v_0.
  Suppose ≥ v_i-1 in the beginning of the current round.
  Because [v_i]=v_i-1 and the non-decreasing property of ,
  we have [+1..v_i]≤ v_i-1≤.
  So [+1..v_i] has been finalized, and we must have >v_i,
  which leads to ≥ v_i in the end of the round.
  Thus the number of rounds is at most k.





To maintain the best decisions, we use an ordered-set of triples B in ([l,r],j) format,
which is the same as the double-ended queue in the sequential algorithm.
However, we still have two steps missing in <ref>:
1) how to find the cordon, and 2) how to update the best decisions.
We now elaborate how to work-efficiently do these two steps.





  §.§.§ Finding the Cordon


To find the ready states in each round, we use the .
Namely, with all states up to  finalized, 
we can attempt to use the tentative states after  to update other tentative DP values. 
Once we find any j that can update i,
we put a  at i. Among all s, the smallest (leftmost) one will give the final position of the cordon. 

However, note that we cannot afford exhaustive checking for all pairs of states (j,i). 
First of all, checking all possible j> may incur large overhead in work,
since most of the later states are unready anyway. 
Ideally, the algorithm should check up to exactly the position of  (but this would be a chicken-and-egg problem). 
To handle this, our idea is to use prefix-doubling (see function FindCordon in <ref>), 
which attempts to extend the cordon by a batch of 2^t-1 states for increasing t in each substep t. 
If the entire batch is ready—i.e., no states in [+1,+2^t) can be relaxed by each other, and all sentinels are outside the batch—we 
try a larger step and extend the cordon to + 2^t+1. 
During the process, we will maintain  as the leftmost sentinel so far. 
Once we find  is inside the batch, 
it means that this batch is not fully ready. 
Therefore, the process stops and returns the current value of  to the main algorithm. 



Using prefix doubling, the parallel algorithm may check more states than the ready ones, 
but the number of “wasted” states is at most twice of the 
“useful” ones which will be finalized in this round. 

Hence, the total number of processed states is O(n). 

We then discuss the way to avoid checking all states i>j when j puts sentinels. 
By DM, if j can successfully relax i, then j can also successfully relax all states i..n.
Therefore, we only need to put a sentinel at the first such state i. 
Recall that we maintain all best decision triples in a data structure B in sorted order.
By DM, we can simply binary search (O(log n) cost) in B to find s_j as the first tentative state that can be updated by j,
and put a sentinel there.


The FindCordon in <ref> gives the full process as described above.
Each iteration of the while-loop at <ref> is a substep,
which processes a batch of states in [+2^t-1,+2^t) in substep t. 
Then for each state j in this batch (in parallel), 
we use B to find the first state that can be updated by j and put a sentinel at this position s_j.
Finally, the leftmost sentinel so far forms the cordon. When the cordon is within the current batch, the algorithm returns. 
We also show an illustration of this process in <ref>. 



  The function FindCordon has O(log n) work and O(log^2 n) span, where =--1 is the frontier size.


As discussed above, the prefix doubling scheme 
  may attempt to process up to ' states, where ' ≤ 2. 
  For each such state, we may binary search in B to find s_j in O(log n) cost, and check the condition on <ref> in O(1) cost.
  Therefore, FindCordon has work O(log n) and span O(log^2 n).





In the , if tentative state j can successfully relax tentative state i,
we should put a sentinel on i.
Then the cordon will be the minimum position of sentinels.
By decision monotonicity, if j can successfully relax i,
then j can also successfully relax any state in [i..n].
So for each j, we can binary search the first position s_j that j can relax s_j,
and put a sentinel on s_j.
As we have stored all best decisions,
we can directly do this binary search on array B,
which takes O(log n) time for each j.



To achieve work efficiency, in each round we need to keep the work of finding the cordon to be proportional to the number of ready states.
So we cannot afford calling the binary search for each j∈[+1,n].
FindCordon shows a dual binary search method.
We keeps doubling the set of ready states
until the minimum sentinel is within the states we have checked,
so all future sentinels will not be smaller than the current one.


Let L=--1,
then FindCordon has O(Llog n) work and O(log^2n) span.





  §.§.§ Generating New Best-Decision Array



In the sequential algorithm, after we compute D[i],
we binary search the first index s_i that i can successfully relax D[s_i..n],
and update [s_i..n] to i.
However, in <ref> we have a batch of newly finalized states D[+1..-1],
and their ranges to update can affect each other.


The efficiency of the algorithm relies on maintaining an ordered data structure B
for all best decision triples. 
We will store B as an array of all such triples in sorted order, such that the binary search 
in <ref> can be performed efficiently. 
Therefore, after we get the newly finalized states D[+1..-1],
we need to update B accordingly to get the new best decision for all states in [, n]. 

We use a divide-and-conquer approach to do this.


Function FindIntervals(,,,) finds all best decision triples for states in range [,], with best decisions
in range [,]. 
Note that we only need the best decisions for all states after . 
All these states must have their current best decisions within [+1, -1] 
(if their best decisions are before , they must have been ready in this round and been included in the frontier). 
Therefore, at the root level, we call FindIntervals(+1,-1, , n). 

In FindIntervals, we first compute =[] where =(+)/2, i.e., the best decision of the state in the middle.
By (convex) DM, the best decisions of i∈[,-1] are in [,], 
and the best decisions of i∈[+1,] are in [,].
We will deal with the two subproblems in parallel.
To collect all ([l,r],j) triples in parallel, we build a tree-based structure bottom-up in the recursion.


Finally, we flatten the tree to an array and 
merge the adjacent intervals if they have the same value of j. 



  The function UpdateBest has O(log n) work and O(log^2 n) span,
  where =--1 is the frontier size. 


  Flattening and removing duplicates can be performed by simple parallel primitives on trees and arrays in O() work and O(log n) span. 
  Below we will focus on the more complicated FindIntervals function. 
  The span of FindIntervals comes from 1) O(log n) levels recursions and 2) O(log n) span to check all states in [,] in parallel.
  For the work, each recursive call in FindIntervals deals with a range of states [,] using best decision candidates in range [,].
  The algorithm first finds ∈ [,] and its best decision ∈ [,]. 
  This can be done by comparing all possible decisions in [,], which is O(-) work. 
  split the ranges into two subproblems and recurse. Let N=|-+1| and M=|-+1| denoting
  the sizes of the two ranges. The work of FindIntervals indicates the following recurrence:
  
    W(N,M)=W(N/2,M_1)+W(N/2,M_2)+O(M), M_1+M_2=O(M)

  This solves to O(Mlog n). On the root level, M=--1=. 
  This proves the stated work bound. 




Our algorithm tackles this by doing the binary search for all ready states at the same time.
Suppose we want to use D[j_l..j_r] to update [i_l..i_r].
We first compute j_0=[i_m].
Then by decision monotonicity, the best decisions of i∈[i_l,i_m-1] is in [j_l,j_0],
and the best decisions of i∈[i_m+1,i_r] is in [j_0,j_r].
After we find all the intervals ([l,r],j),
we merge the adjancent intervals if they have the same j.



It can be shown that the work of UpdateBestChoice is proportional to the number of ready states.
The only problem is how to efficiently adding new intervals into B while keeping their order.
As the function FindBestChoiceIntervals is called recursively in two directions,
we can build a binary tree of intervals based on the recursion structure,
and then flatten the tree into an array.
This method does not bring any race condition.
Another method is to maintain the min/max index for each ready state j.
When we add a new interval ([l,r],j), we call one write_min and one write_max for state j.
At last we can easily build B based on the range of each j.

Let L=--1,
then <ref> has O(Llog n) work and O(log^2n) span.




 §.§ Parallel Concave GLWS





To extend the algorithm to the concave case, we need a few modifications.
In FindCordon, by the concavity,
if j can update i, then j must be able to update j+1.
Therefore, in <ref> in <ref>, we check whether j can update j+1. 
If so, we put a  at j+1.
The other modifications are in FindIntervals. 
First, due to concavity,  
when we find  as the best decision of  in <ref>,
we need to swap the last two parameters in the first and second recursive calls, 


i.e., the best decision range for states +1 to  must be before , and those in  to -1 must be after . 

A more involved modification in the concave GLWS is that after we get the array B from FindIntervals,
we have to merge it with the old array B before this round—FindIntervals only considers the best decisions among [+1,-1],
but in the concave case, these states may also have better decisions using states before .
Suppose we have generated the array B_ storing ([l,r],j) triples,
and we want to merge it with B_.
Both of B_ and B_ contain the best decisions of states [..n].
The difference is that the js in B_ are from [0..],
while the js in B_ are from [+1..-1].
By the concave decision monotonicity,
the key is to find a cutting point p, where
the best decisions of [..p] are from B_,
and the best decisions of [p+1,n] are from B_.



Here we show <ref> to find p in O(hlog n) work and O(log n) span,
where h=|B_| is the frontier size.
For all l_k in B_, we pre-process its best decision stored in B_.
This step requires O(hlog n) work and O(log n) span.
Then we search in B_ to find the interval that p locates in.
After this step, there is only one interval ([l_k,r_k],j_k) in B_ is interesting.
Then we can binary search in B_ to find the exact p.
Note that this method can be easily modified to merge B_ and B_ even if the cost function is convex.
For page limitation, we elaborate on this part in <cit.>.
With careful design, this part can also be finished in O(log n)
work and O(log n) span, where  is the frontier size.



However, this is not a perfect parallelization.
The number of rounds may be larger than the depth of the perfect DAG.




 §.§ Theoretical Analysis


In this section, we show theoretical analysis for our parallel GLWS algorithm. We first summarize our main results as follows. 



  
  Given an input sequence of size n, the  for the convex GLWS has O(nlog n) work and O(klog^2n) span,
  where k is the effective depth of the perfect DAG 
  G^*_ of algorithm  introduced in <ref>. 
  It is a perfect parallelization of . 
  
  
  Given an input sequence of size n, and the sequential GLWS algorithm 
   introduced in <ref>,   
  the  for the convex GLWS has O(nlog n) work and O(klog^2n) span,
  where k is the effective depth of the perfect DAG 
  G^*_. 
  It is a perfect parallelization of . 
  
  Given an input sequence of size n, and the sequential GLWS algorithm 
   introduced in <ref>,   
  let k=(G^*_) be the effective depth of the 
  -perfect DAG. 
  Then the  for the convex GLWS has O(nlog n) work and O(klog^2n) span. 
  It is a perfect parallelization of .  



More intuitively, k in <ref> is also the number of best decisions to make in the final solution: 
for the post-office problem, it is the number of post offices in the optimal solution. 




  The  for the concave GLWS has O(nlog n) work and O(klog^2n) span,
  where n is the input size and k is the effective depth of the optimized DAG G_
  algorithm introduced in <ref>. 
  It is an optimal parallelization of .   
  
  Given an input sequence of size n, and the sequential GLWS algorithm 
   introduced in <ref>,
  the  for the concave GLWS has O(nlog n) work and O(klog^2n) span,
  where k is the effective depth of the optimized DAG G_
  algorithm introduced in <ref>. 
  It is an optimal parallelization of .    
  
  Given an input sequence of size n, and the sequential GLWS algorithm 
   introduced in <ref>,   
  let k=(G_) be the effective depth of 
  
  the -optimized DAG. 
  Then the  for the concave GLWS has O(nlog n) work and O(klog^2n) span. 
  It is an optimal parallelization of .  



We first prove that both algorithms are nearly work-efficient and have O(nlog n) work.



  The  for GLWS has O(nlog n) work for both convex and concave case.


  Combining <ref> (and the discussion in <ref>), the work for each round is O(log n), 
  where =--1 is the frontier size. 
  Since the frontier sizes  across all rounds add up to n, 
  the entire algorithm has O(nlog n) work. 
  
  


We then show that the number of rounds in both convex and concave cases is the effective depth of G_.
Recall that the DAG G_ includes normal edges between all states j and i<j,
and effective edges between a state j and its best decision. 
The effective depth (G_) is the largest number of effective edges in any path. 


  The  for GLWS finishes in k rounds, where k is (G_). 


  Define the effective depth (s) of a state s
  as the largest number of effective edges of a path ending at s. 
  We will inductively prove that a state s is in the frontier of round r iff. s has effective depth r.
  The base case (boundary cases) holds trivially.
  
  
  Assume the conclusion is true for r-1. We first prove the “if” direction, i.e., if a state s has effective depth r,
  it must be in the frontier of round r. This is equivalent to show that there is no sentinel on all states from  to s.
  Assume to the contrary that there is a state y∈ (,s] with a sentinel, which is put by state x∈ (, y].
  This means that x is a better decision for y than all states before , 
  indicating that y's best decision y^*≥.
  
  Based on the induction hypothesis, the effective depth of y^* must be larger than r-1. 
  Therefore, (y)= (y^*)+1 >r-1+1=r, which means that (y) is at least r+1.
  Based on the recurrence, there is a normal edge from y to s, so (s)≥ r+1, leading to a contradiction.
  
  We then prove the “only if” condition, i.e., if a state s is in the frontier of round r, it must have effective depth r. 
  The induction hypothesis suggests that all states with effective depth smaller than r have been finalized in previous rounds, so we only
  need to show that (s) cannot be larger than r. 
  Assume to the contrary that (s) ≥ r+1. 
  Let the path to s with effective depth (s) be x_1, x_2, …, s. 
  Since the total number of effective edges on this path is at least r+1, there must exist
  an effective edge x_ix_i+1 on the path such that (x_i)=r and (x_i+1)=r+1.
  However, based on the induction hypothesis, x_i's best decision must have been finalized. 
  During <ref>, x must get its true DP value, and will find itself able to update x_i+1.
  Therefore, there will be a sentinel on x_i+1≤ s, and s cannot be identified in the frontier of round r.


We will then show that the number of rounds of the convex case is also the effective depth of the 
-perfect DAG G^*_. 
This is stronger than the -optimized DAG as shown above. 
Recall that the perfect DAG G^*_ contains all best decision edges in G_. 



  
  
  The  for convex GLWS runs in k^* rounds, where k^* is (G^*_). 


  Define the perfect depth (s) of a state s
  as the largest number of effective edges of any path ending at s in G^*_.   
  Similarly, we will show by induction that in round r, all states with perfect depth r will be processed. The base case holds trivially. 
  Assume the conclusion holds for r-1. 
  In round r, we will show that a state s with perfect depth r must be put in the frontier. 
  Let s^* be the best decision of s, then (s^*)=r-1 and therefore s^*<. According to DM, any state x between  and s must have its best decisions
  x^*≤ s^*<, indicating that (x^*)≤ r-1. Therefore, x must find its true best decision in B, and cannot be updated by any other tentative states in <ref>. This means that there will be no sentinel between  and s, so s must be identified ready in round r. 
  Therefore, a state with perfect depth r must be finalized in round r, leading to the stated theorem. 



Combining <ref>

proves the span bounds in <ref>.


Discussions with totally monotonicity.
A more general condition to enable DM is convex/concave totally monotonicity <cit.>.
Our algorithm works as long as the transition function is convex/concave totally monotone. 
For page limit, we formally explain this in <ref>full version of this paper <cit.>.






  The  for the convex GLWS has O(nlog n) work and O(klog^2n) span,
  where n is the problem size and k is the longest dependency path. It is a perfect parallelization of the Γ_1D
  algorithm introduced in <ref> for the convex case. 






§ OTHER PARALLEL DP ALGORITHMS


We now show that our algorithmic framework can be used to parallelize a wide variety of classic sequential DP algorithms.

Due to the space limit, we will focus on two of them and leave others in the appendicesthe full version of this paper <cit.>, with a short summary in <ref>.


In particular, for the optimal alphabetic tree (OAT) problem (<ref>), we partially answered a long-standing open problem by <cit.> for reasonable input instances (for instance, positive integer weights in range n^(n)).
For the GAP problem (<ref>), we showed the first nearly work-efficient algorithm with non-trivial parallelism.
More interestingly, this algorithm combines all techniques in the algorithms for convex GLWS and sparse LCS. 







 §.§ Parallel Optimal Alphabetic Trees (OAT)


The optimal alphabetic tree (OAT) problem is a classic problem
and has been widely studied both sequentially <cit.>
and in parallel
<cit.>.
Given a sequence of non-negative weights a_1..n, the OAT is a binary search tree with n leaves and has the minimum cost,
where the cost of a tree T is defined as:

    cost(T)=∑_i=1^na_id_i

Here d_i is the depth of the i-th leaf (from the left) of T (the root has depth 0).
One can view the weight a_i as the frequency of accessing leaf i, and the depth of a leaf is the cost of accessing it.
Then the cost of T is the total expected cost of accessing all leaves in T. 


The OAT problem is closely related to other important problems such as the optimal binary search tree (OBST) <cit.>
and Huffman tree <cit.>. 



Sequentially, <cit.> showed an OAT algorithm with O(nlog n) work. 
Later <cit.> simplified this algorithm. 
In parallel, <cit.> showed an algorithm based on Garsia-Wachs.
We will apply our techniques to this algorithm to improve the span bounds. 
Due to page limit, we provide the details of  <cit.> in <ref>the full version paper <cit.>, and review the high-level idea here.
The algorithm computes an l-tree <cit.>,
which has the same depth with and will be finally converted to the OAT in O(n) work and polylogarithmic span. 
The key insight of  <cit.> is to start with a sequence of n leaf nodes with the input weights, 
and find several disjoint intervals in the sequence to process in parallel. 
This partition is done by various operations on the Cartesian tree of the input sequence, which requires O(nlog n) work and O(log^2 n) span.
Larmore et al. showed that processing each interval can be reduced to a convex LWS. 
The solution of the LWS will connect items in this interval into a forest, 



which becomes a subgraph in the final l-tree. 








Finally, for each tree in the forest, we insert its root back to the sequence and repeat the process. 
This reinsertion step takes O(nlog n) work and O(log n) span by basic parallel primitives such as sorting and range-minimum queries. 

The further rounds will connect the forest to the final l-tree. 



Larmore et al. also showed that the number of such intervals shrinks to half in each iteration, so the algorithm will finish in O(log n) rounds. 
Here all other steps in addition to solving convex LWS take O(nlog^2 n) work and O(log^3 n) span. 
The remaining cost of the algorithm is to solve convex LWS in each round, multiplied by the number of rounds, which is O(log n).

Larmore et al. originally used the parallel convex LWS algorithm from <cit.>, 
which has O(m^2log m) work and O(log^2m) span when taking an input interval of length m. 


Later, <cit.> improved the parallel convex LWS algorithm to O(m^1.5log m) work and O(√(m)log m) span,
yielding O(n^1.5log^2n) work and O(√(n)log^2n) span for the OAT algorithm—the work overhead is still polynomial.



<cit.> left the open problem on whether there exists an OAT algorithm with Õ(n) work and (n) span, 
which remains unsolved for three decades.

Note that the convex LWS problem is a special case of the convex GLWS problem discussed in <ref> with E[i]=D[i].
Hence, <ref> directly gives O(mlog m) work and O(klog^2m) span for convex LWS problem, and here k is the longest dependency path of best decisions.




In Larmore's algorithm, the forest for each interval is constructed iteratively by the DP algorithm on LWS: 
if iteration i finds the best decision at iteration j<i, 
then iteration i creates one more level on top of the forest at iteration j.
This means that the k is equivalent to the depth of the forest, which is upper bounded by the final OAT height h. 
We present more details in <ref><cit.>.


Hence, we can parameterize our final bounds using h as:



  The optimal alphabetic tree (OAT) can be constructed in O(nlog^2n) work and O(hlog^3n) span,
  where n is the size of input weight sequence and h is the height of the OAT.



This algorithm is nearly work-efficient with span parameterized on h.
One useful observation is that the OAT height h is polylogarithmic with real-world input instance with positive integer weights and fixed word length.

More formally, we can show that:



  In an OAT, the subtree weight grows by at least twice for every three levels.





  If all input weights are positive integers in word size W, 
  the OAT height is O(log W). 



The proof is not complicated and we provide it in <ref><cit.>.
With this lemma, we can state the following corollary.



  If the input key weights are positive integers with word size W=n^ n,
  the OAT can be constructed in O(nlog^2n) work and (n) span,
  where n is the input size.





The bounds also hold for real number weights if the ratio between the largest and smallest weight is n^ n.
We note that in realistic models we usually assume word-size W=n^O(1), 

in which case <ref> affirmatively answers the open problem in <cit.>.

 



 §.§ The GAP Edit Distance Problem







The GAP problem is a variant of the famous edit distance problem.
The GAP problem aligns two input strings with sizes n and m≤ n, and allows editing a substring with certain cost function (see formal definition below).
This problem has been widely studied both sequentially <cit.>
and in parallel <cit.>.

As noted by <cit.>, most real-world cost functions are either convex or concave, yielding Õ(nm) work for the GAP problem sequentially.
Unfortunately, to the best of our knowledge, these existing parallel algorithms for the GAP problem need Ω(n^2m) work,
and the polynomial overhead makes them less practical.

More specifically, GAP takes two strings A[1..n] and B[1..m], and computes
the minimum cost to align A and B using the following operations:
1) deleting A[l+1..r] with cost w_1(l,r),
and 2) deleting B[l+1..r] with cost w_2(l,r).
Here we consider the following recurrence, which is usually referred to as the GAP recurrence:

    P[i,j]   =min_0≤ i'<iD[i',j]+w_1(i',i) 
    
    Q[i,j]   =min_0≤ j'<jD[i,j']+w_2(j',j) 
    
    D[i,j]   =min{P[i,j],Q[i,j], D[i-1,j-1] | A[i]=B[j]}.

Here P[i,j] and Q[i,j] indicate the edits on the two strings.
Directly computing the recurrence uses O(n^2m) work.
Since most real-world cost functions in machine learning, NLP, and bioinfomatics <cit.> are either convex or concave, 
sequentially each row in P or column in Q is a convex or concave GLWS and can be computed in Õ(n) or Õ(m) work.
Hence, computing the entire P and Q takes Õ(nm) work, leading to the same cost for computing D and the entire problem. 
We denote this standard sequential algorithm as . 



Parallelizing this approach is extremely challenging even with the parallel convex/concave GLWS in <ref> as a subroutine, and we are unaware of any existing work on this.
The challenge here is that the rows in P interact with the columns in Q. 
For instance, computing a row in P requires one element from each column in Q, but computing those elements again requires previous rows in P.

Our key insight to parallelize this algorithm is to use the  to efficiently mark the ready region to be computed in each round.
Note that as a generalization of the classic edit distance/LCS, the GAP recurrence is similar to Recurrence <ref>, but with “jumps” in computing P and Q.
An illustration is given in <ref>.
In addition to the diagonal edges as in LCS (see <ref>), for rows and columns, 
there also exist effective (red) edges (see <ref>(a)). 
Here for simplicity we only draw a subset of these edges, 
and every state D[i,j] always have one vertical effective edge (to compute Q[i,j]), one horizontal effective edge (to compute P[i,j]), and may have a diagonal edge if A[i]=B[j].
All these edges imply the sentinels, which form the cordon and imply
the regions for ready states, as shown in <ref>(b). 
The cordon is still a staircase as in LCS. 

However, finding the cordon in GAP is sophisticated.
We cannot directly use a tournament tree as in LIS, since the vertical and horizontal edges are computed on-the-fly and not known ahead of time.
Meanwhile, in a 2D table where the cordon is a staircase, we cannot simply use prefix-doubling as in GLWS in <ref>.
We propose a unique solution here to use prefix-doubling on a 2D table and computes the staircase cordon efficiently.
This approach will consider each row separately, but for all rows, we run prefix doubling synchronously and try to see if the next ranges are available.
First, we put a  on state (x,y) with a diagonal edge if (x-1,y-1) is not finalized. 
We will maintain the best-decision structure for each row and column, in the same way as the GLWS algorithm. 
For this region to be checked, we will use the same approach as in <ref> to compute P and Q, take the minimum as D, and use D to check their readiness.
If a state (x,y) obtains the best decision from another tentative state, we will put a  on (x,y), 
which will block the other states (x',y') with x'≥ x and y'≥ y.

The work to put the s is proportional to the number of states we checked in the prefix-doubling, and the span is polylogarithmic.

Finally we discuss how to handle the s placed as above.
We store all s based on the row index on increasing order.
After this, applying a prefix-min on these s gives part of the cordon (if they exist), and we will merge it with the previous cordon.
Then, for all tentative states, we check whether they are on the correct side, and invalidate those across the cordon. 
Since we are using prefix doubling, the wasted work for the invalid states can be amortized.
In the next prefix doubling step, we will also use the cordon to limit the search region.
Once all states within the cordon are checked for readiness, we can move to the next round.
Due to prefix doubling, we only need O(log n) steps in each round.



  The  for the GAP problem has O(mnlog n) work and O(klog^2n) span,
  where n and m≤ n are the input size and k is the effective depth of the -optimized DAG
  for the sequential algorithm  introduced in <ref>.





For page limitation, we present the proof in <cit.>.
The work bound can be inferred from the O(nlog n) work for each 1D GLWS problem.
The proof of span bound is similar to <ref>.




Proof of <ref>.
Recall that the sequential GAP algorithm  gets the DP value
for each state s=(i,j) by solving the GLWS problems in row i
and column j, respectively, and the diagonal edge (i-1,j-1)(i,j) if applicable.
Therefore, the optimal DAG G_ contains three types of edges 

  
  * (i,j)(i',j) for all i'>i,
  
  * (i,j)(i,j') for all j'>j, and
  
  * (i-1,j-1)(i,j) if A[i]=B[j]. 

Among them, the effective edges include:

  
  * (i,j)(i',j) where i is the best decision for i' in the GLWS problem on row j,
  
  * (i,j)(i,j') where j is the best decision for j' in the GLWS problem on column i, and
  
  * (i-1,j-1)(i,j) if A[i]=B[j]. 
 

WLOG we assume m≤ n in this section. We first prove the span bound. We will show that the  finishes in k rounds, where k is the effective depth of G_. 

  Given two sequences of sizes n and m≤ n,
  the  on GAP edit distance finishes in k=(G_) rounds. 


  The proof is similar to <ref>. We also define the effective depth of a state s as s.
  We will show by induction that s is processed in round r iff. s=r. The base case holds trivially. 
  
  Assume the conclusion holds for all rounds up to r-1. We will show it is also true for round r. 
  We first prove the “if” direction, i.e., if a state s=(i,j) (i-th row, j-th column) has effective depth r,
  it must be in the frontier of round r. 
  This is equivalent to show that there is no sentinel that blocks s.
  For simple description, for two states s=(i,j) and s'=(i',j'),
  we say s≺ s' if i≤ i' and j≤ j'. 
  Clearly, if a state s≺ s', a sentinel on s will block s'. 
  Assume to the contrary that there is a state y≺ s with a sentinel, which is put by another tentative state x≺ y.
  This means that the tentative state x is a better decision than all finalized states, indicating that the best decision of y, denoted as y^*, must 
  also be tentative.
  
  Based on the induction hypothesis, the effective depth of y^* must be larger than r-1. 
  Therefore, (y)≥(y^*)+1 >r-1+1=r, which means that (y) is at least r+1.
  
  Let y=(i',j'). y and s can be connected by either one normal edge (when they are in the same row or column) 
  or two normal edges ((i',j')(i,j') and (i,j')(i,j)). 
  This means that the effective depth of s is at least the same as y, which is r+1. 
  This leads to a contradiction. 
  
  We then prove the “only if” condition, i.e., if a state s is in the frontier of round r, it must have effective depth r. 
  The induction hypothesis suggests that all states with effective depth smaller than r have been finalized in previous rounds, so we only
  need to show that (s) cannot be larger than r. 
  Assume to the contrary that (s) ≥ r+1. 
  Let the path to s with effective depth (s) be x_1, x_2, …, s. 
  Since the total number of effective edges on this path is at least r+1, there must exist
  an effective edge x_ix_i+1 on the path such that (x_i)=r and (x_i+1)=r+1.
  However, based on induction hypothesis, x_i's best decision must have been finalized. 
  In round r, x must get its true DP value, and will find itself able to update x_i+1.
  Therefore, there will be a sentinel on x_i+1≺ s, and s cannot be identified in the frontier of round r.


Combining <ref>, the span in each round is O(log^2n).
This proves the span bound in <ref>. 

We then prove the work bound in <ref>. 

  Given two sequences of sizes n and m≤ n,
  the  on GAP edit distance has work O(mnlog n). 


  As G_ is a grid graph, its depth is no more than m+n.
  By <ref> the algorithm will finish in k=O(n) rounds.
  In each round, we do prefix-doubling across all m rows and try to push the frontier on each row.
  In each prefix-doubling step we do a prefix-min that costs O(m) work,
  so the cost of prefix-doubling is O(mlog n) in each round, and O(mnlog n) in total.
  Suppose h is the frontier size in one round.
  Due to prefix-doubling, the number of tentative states we visited is at most 2h.
  Combining <ref>,
  in each row/column we can achieve work proportional to the number of tentative states.
  Thus the cost to put sentinels and maintain the best decision arrays is also O(mnlog n).






The computation of D[i,j] reduces to the computation of P[i,j] and Q[i,j].
In each column/row, the computation of P[i,1..n] and Q[1..m,j] is a GLWS problem.
The step to take the minimum of P[i,j], Q[i,j] and D[i-1,j-1] (if possible)
plays the role as the f() function in GLWS.
If the weight function w_1 and w_2 is convex/concave,
a sequential algorithm can make use of GLWS algorithm with DM to compute all D[i,j] in O(mnlog mn) time.

To parallelize the 2D GAP algorithm,
we can borrow the staircase idea from the parallel LCS algorithm in <ref>.
The difference is that in addition to the matching positions (A[i]=B[j]),

the sentinels are put based on the GLWS problems of row i and column j.
By  we need to put a  on all state that relies on any previous tentative state.
An example is shown in <ref>,
where D[3,4] can relax P[3,7], and D[3,3] can update Q[5,3].
So we put two s at (3,7) and (5,3), plus the  at (6,1) which is a matching position.
Every  will block the upper-right area of itself,
so we can build the cordon as <ref>.


If w_1 and w_2 is convex/concave, we can achieve work-efficiency by the prefix-doubling technique.
In each round of , we try to push the frontier of each row in 1,2,4,8 steps.
We can check the inner dependencies with work proportional to the number of tentative states we have visited.
During this process, if we put a  at (x,y),
we set a limit that the frontier in these rows cannot be pushed beyond y.

The analysis of the span of our algorithm is as follows.
Let _P[i,j] and _Q[i,j] be the best decision of P[i,j] and Q[i,j], respectively.
In the begining the Γ-optimized DAG is the grid.
For each state (i,j),
we add two effective edges: from _P[i,j] to (i,j), and from _Q[i,j] to (i,j).
And if A[i]=B[j], we add an effective edge from (i-1,j-1) to (i,j).
Then the number of rounds of our algorithm is the longest path of effective edges.







 §.§ Other Algorithms







We further apply the  to several other problems in <ref>the full version of this paper <cit.>, including 
a GLWS problem on a tree structure, 
k-GLWS, and Optimal Binary Search Tree (OBST). 
For both k-GLWS and OBST, our algorithm leads to a correct algorithm with trivial parallelism. 
However, we show that they are still optimal parallelization of the sequential algorithms—this indicates that to further improve the parallelism, 
one should probably find or redesign another sequential algorithm to start with.







 §.§ General LWS on Trees





The idea of decision monotonicity (DM) can be applied to various structures more than just 1D cases discussed in <ref>.
The efficient parallelism on 2D grid structure is introduced in <ref>, and we now show the techniques to enable high parallelism on the tree structure.
Here we refer to this problem as .


Let T be a tree with n+1 nodes, and node 0 is the root.
We use p(v) to denote the parent of node v and d_v be the distance from node 0 to node v.
 takes the input tree T, a cost function w, and the boundary D[0], and computes:

    D[v]=min{E[u]+w(d_u,d_v)}

where u is any ancestor of v, and E[u]=f(D[u],u) that can be computed in constant time from D[u] and u.
The cost function w is decided by the depths of u and v.
Note that here sibling nodes v_1 and v_2 will have the same DP value,
but E[v_1] and E[v_2] can be different given that v_1 and v_2 are also part of the parameter in computing the function f.

In this section we assume w is convex, but our algorithm can adapt to the concave case with some modifications.



  §.§.§ Building Blocks

We will first overview some basic building blocks, which are crucial subroutines used in our algorithm.


Persistent Data Structures.
A persistent data structure <cit.> keeps history versions when being modified.
We can achieve persistence for binary search trees (BSTs) efficiently by path-copying <cit.>, where only the affected path related to the update is copied.
Hence, the BST operations can achieve persistence with the same asymptotical work and span bounds as the mutable counterpart.

Heavy-Light Decomposition (HLD).
HLD <cit.> is a technique to decompose a rooted tree into a set of disjoint chains. 
In HLD, each non-leaf node selects one heavy edge,
the edge to the child that has the largest number of nodes in its subtree.
Any non-heavy edge is a light edge.
If we drop all light edges, the tree is decomposed into a set of top-down chains with heavy edges.
As such, HLD guarantees that the path from the root to any node v contains O(log n) distinct chains plus O(log n) light edges.
If we use BSTs to maintain each heavy chain in HLD, we can answer path queries (e.g., query the minimum weighted node on a tree path) in O(log^2n) work.

Range Report Based on Tree Depth.
We now discuss a data structure that efficiently reports the set of nodes in a subtree of T where the depths of the nodes are in a given range l to r.
First we build the Euler-tour (ET) sequence of T, so any subtree of T will be a consecutive subsequence in the ET.
We can map all nodes to a 2D plane each with coordinates (f_v,d_v), where f_v is the first index of v in the ET, and the d_v is the tree depth of v.
Now the original query is a 2D range report on this 2D plan.
A range tree <cit.> can be built in O(nlog n) work and O(log^2 n) span, and answer this query in O(m+log^2n) work and O(log^2n) span where m is the output size.



  §.§.§ Our Main Algorithm


Here if we consider any tree path, Recurrence (<ref>) is exactly the same as for the 1D case in <ref>.
Hence, we can use a similar approach as in <ref> by maintaining the best-decision array of ([l,r],j) triples,
meaning that for elements v with depth l≤ d_v≤ r, v's best decision is j.
The challenge here is the branching nature of a tree—we need to handle path divergences at nodes with more than one child.
The work can degenerate to Õ(n^2) if we copy the best-decision arrays at the divergences, since we can end up with O(n) leaves.
Sequentially, we can depth-first traverse the tree and compute the “current” best-decision array, and we only need to revert the array when backtracking.
However, this approach is inherently sequential.
To utilize the  on a tree structure, we need to resolve the following two challenges: 1) how to efficiently identify the ready nodes;
and 2) how to efficiently maintain the best decision arrays for each node.




Identifying the Ready States.
Similar to the 1D case <ref>, we maintain the best-decision array for each tree path. 
Then we traverse the tree top-down, and we identify the ready states v that can be finalized in this round, and compute them in parallel.
An illustration can be found in <ref>.


Our high-level idea still follows the prefix-doubling technique, similar to the 1D case.
In the t-th doubling step we expand all nodes with 2^t-1≤ d_v < 2^t.
These nodes can be extracted by a range report shown in <ref>. 
We use prefix doubling and the checking process in <ref> to decide the boundary that forms the cordon in the next round.

When checking the availability, we can use the HLD-based tree path query to find the minimum (highest) node on each path that is not available.
We will put sentinels on these nodes that block their subtrees.
The process stops when we find such nodes for all tree paths. 
In <ref>, these ready nodes are shown in green.
In the next round we can asynchronously work on the subtrees on the cordon in parallel.
We repeat this process until all nodes are finalized (correctly computed).

Here one difference to the 1D case is that the work cost of perfix doubling cannot be perfectly amortized.
In the 1D GLWS, if the prefix-doubling stops at step t,
we visit at least 2^t-1 ready states and at most 2^t-1 unready states,
thus the work to visit the unready states can be amortized.
However, in the tree case we are doing prefix-doubling by the depth of nodes.
The number of nodes in the last prefix-doubling step can be much larger than in the previous steps,
and the cost cannot be amortized.
The insight is that due to the prefix-doubling,
each node v will be visited in at most O(log n) rounds.
Plus the O(log^2n) work of the range report,
the work in each round can be amortized to O(hlog^3n), where h is the frontier size.


Updating the Best-Decision Arrays.
The most interesting part in this algorithm is how to maintain the best-decision arrays for all tree paths while achieving work efficiency and high parallelism.

Due to the tree structure, the best-decision arrays for different branches of the tree share some parts.
In total, there can be O(n) paths with total sizes of O(n^2).
The key challenge is to save the work and space by sharing parts of the arrays, while updating them highly in parallel.


Consider the simple case when the ready nodes form a chain (the 1D case in <ref>).
Here we use persistent BSTs to maintain the best-decision arrays on each node.
We first use UpdateBestChoice in <ref> to generate the best-decision array B in the middle node of the chain,
and merge it with the old B (the one stored at the node above this chain) using the similar technique in <ref>.
During this process we use path-copying to generate a new version the new array.
Then we work on the upper part and the lower part of the chain in parallel.
By this divide-and-conquer method, we can generate the best-decision array on each node of the chain
with O(mlog^2 m) work and O(log^3m) span, where m is the length of the chain.

In the general case, the structure of ready states can be arbitrary.
To achieve work-efficiency and high parallelism, our solution is in a “BFS-style” algorithm that utilizes the properties of HLD (see <ref>). 
For all ready nodes in this round, we extend the heavy chain that is directly connected to the finalized nodes.
Since the heavy chain will not diverge, the approach is the same as the 1D case except for additional persistence, with work proportional to the total number of nodes and polylogarithmic span.
Once we finish updating the heavy chain, we will in parallel work on the light children of the nodes on the heavy chain we just proceeded.
The overall structure is similar to a BFS with heavy edges with weight 0 and light edges with weight 1.
Since each node only appears in one heavy chain, the work is still proportional to the number of ready nodes.

We can also achieve high parallelism due to the fact that there can be at most O(log n) heavy chains and light edges from the root to any node v.
Hence, we can finish updating all paths in a logarithmic number of steps per round, which guarantees both work-efficiency and high parallelism.

Here we assume we build the HLD for the entire tree T at the beginning, but we can also build the HLD for the ready states locally for each round.

Combining all pieces together, 
in each round we can determine the ready states and maintain the best-decision arrays
with work proportional to the number of ready states and polylogarithmic span.
We hence have the following theorem:

   solves  in O(nlog^3n) work and O(klog^4n) span,
  where k is the longest path in the best decision dependency graph.
 


For convex GLWS on tree, the best decisions of nodes in a subtree is still monotone with the depth, just like the convex 1D GLWS.
So for each subtree we can still maintain the best decision array of ([l,r],j) triples,
meaning that in this subtree if l≤ d_v≤ r then v's best decision is j.
However, due to the branches in the tree, the best decision array needs to be persistent.
Each node v store a version of the best decision array,
which is the same as the best decision array if we extract the path from 0 to v as a prefix of the 1D problem.
Sequentially, we can do a breath-first-search (BFS) or depth-first-search (DFS) on the tree while maintaining the best decision array on each node.
In addition, if we do DFS, the best decision array does not need to be persistent, but only needs to be revertable,
alghouth this method is off our topic in parallelization.

So we consider the simple case in <ref> where state 0 is finalized and state 1 is tentative.
The tentative DP values in this subtree can be computed from the best decision array stored at node 0.
The goal is to determine the maximal connected subgraph that is connected to 1,
in which all states are ready (the green nodes).

For each new node v, we do a binary search to find s_v,
which means v can relax those nodes with depth no less than s_v.
Then v can put a at depth s_v in its subtree.
A node v is ready if for any node u on the path from 1 to v, s_u>d_v, otherwise v is blocked.
This is where we need the tree path minimum query.
The process stops when all nodes in the frontier are blocked.
Then we have a connected tree structure of all ready states
(the green nodes in <ref>).
After we update the best decisions, we can work on deeper subtrees in parallel.


 



 §.§ Parallel k-GLWS


Another well-known variant of GLWS is to limit the output that contains a fixed given number of k clusters in the output <cit.>.
Here we refer to it as the k-GLWS problem.

Formally, let D[i,k'] be the minimum cost for the first i elements in k' clusters, and the DP recurrence is:

    D[i,k']=min_j<iD[j,k'-1]+w(j,i)

where w(j,i) is the cost of forming a cluster containing elements indexed from j+1 to i, and the boundary case D[0,0]=0 and D[i,0]=+∞ for i>0.
Directly solving this recurrence takes O(kn^2) work.
When the cost function w is convex (which happens in many practical settings), the computation of each column in the DP table is a static matrix-searching problem, i.e., for a totally monotone matrix A where A[j,i]=D[j,k'-1]+w(j,i),
we want to compute the minimum element in each column of A.
Theoretically this problem can be solved in O(n) work by the SMAWK algorithm <cit.>, but this algorithm is quite complicated and inherently sequential.
Practically, there exists a simple divide-and-conquer algorithm with O(nlog n) work <cit.>, which is similar to the function FindIntervals in <ref>.
This algorithm first computes the minimum element in the (n/2)-th column by enumerating all elements in this column, and recurse on two sides.
Due to the monotonicity, the minimum element in the (n/2)-th column limits the searches on both sides and guarantees the search ranges shrink by a half.
Hence, the work spent on each recursive level is O(n), yielding O(nlog n) total work for a recursive structure with log n levels.
By parallelizing the divide-and-conquer and using parallel reduce to find the minimum element (with O(log n) span), the total span is O(log^2 n).

We now show that when applying  to this sequential algorithm, the k'-th frontier contains all states D[·,k'].
We can see that in the first round, states D[·,1] are ready.
Since all states D[·,2] depend on some state from D[·,1], we will put s on all states D[·,2] and they thus block all later states.
Then we can inductively show that this applies to all rounds, so finishing this computation requires k rounds.
Then computing all states D[·,k'] in each round using the aforementioned algorithm requires O(nlog n) work and O(log^2 n) span.
Hence, the entire algorithm has O(knlog n) work and O(klog^2n) span.
In this problem, k is also the depth of the DP DAG, so this algorithm is a perfect parallelization of the classic sequential algorithm.


Here we show an parallel algorithm that can be directly derived by the  framework.
In the beginning we finalize all D[i,0] by the boundary condition.
By the , in each round we need to determine the ready states,
and use their DP values to update other tentative states.
As D[*,k+1] relies on D[*,k], we will put sentinels on column k+1,
thus the ready states are actually the column D[*,k].
Then we try to use the D[*,k] to relax D[*,k+1].
We can use the similar divide-and-conquer technique as the FindIntervals function in <ref>
to achieve O(nlog n) work and O(log^2n) span.
The whole algorithm takes O(knlog n) work and O(klog^2n) span.
In this problem k is also the depth of the DP DAG, so this algorithm is a perfect parallelization.
 



 §.§ Optimal Binary Search Tree (OBST)


(Static) OBST is one of the earliest examples of DM optimization.
Given an array of frequency a_0..2n, it computes the recurrence 

    D[i,j]=min_j≤ k<iD[j,k]+D[k+1,i]+w(i,j)

where w(i,j)=∑_k'=2(i-1)^2ja_k', and returns D[1,n].
<cit.> first showed that computing this recurrence only needs O(n^2) work, and later <cit.> showed that this algorithm applies to any convex function w(i,j).
Here let the best decision of a state D[i,j] be the index k that minimizes D[i,j] in <ref>.
In this algorithm, D[i,j] depends on D[i,j-1] (let l be its best decision), D[i+1,j] (let r be its best decision), D[i,l..r], and D[l..r,j].
When applying , due to the dependence from D[i,j] to D[i,j-1] and D[i+1,j], the δ-th frontier contains the states D[i,i+δ].
Hence, although it results in optimal parallelization to the standard sequential algorithm, the algorithm requires n-1 rounds and thus has O(nlog n) span. 
Achieving o(n) span may need new insights to redesign the dependencies. 




GLWS on trees.
We discuss how to apply the  to a GLWS problem on a tree structure in <ref>the full version of this paper <cit.>, which requires a list of algorithmic components such as the heavy-light decomposition, persistent data structures, range trees, and the Euler tour. 
We show that our parallel algorithm is nearly work-efficient (Õ(n) work) and has perfect parallelism. 

k-GLWS.
k-GLWS is another commonly seen problem that can be solved by DP. 
It is similar as GLWS problem as in <ref>, but requires a fixed number of k clusters 
(e.g., exact k post office in the running example). 
In this case, let D[i,k'] be the minimum cost for the first i elements in k' clusters, and the DP recurrence is D[i,k']=min_j<iD[j,k'-1]+w(j,i).
In <ref>the full version of this paper <cit.> we show that applying  here gives an algorithm with Õ(nk) work and Õ(k) span.

Optimal Binary Search Tree (OBST).
(Static) OBST is one of the earliest examples of DM optimization.
Given an array of frequency a_0..2n, it computes the recurrence 

    D[i,j]=min_j≤ k<iD[j,k]+D[k+1,i]+w(i,j)

where w(i,j)=∑_k'=2(i-1)^2ja_k', and returns D[1,n].
<cit.> first showed that computing this recurrence only needs O(n^2) work, and later <cit.> showed that this algorithm applies to any convex function w(i,j).
Here let the best decision of a state D[i,j] be the index k that minimizes D[i,j] in <ref>.
In this algorithm, D[i,j] depends on D[i,j-1] (let l be its best decision), D[i+1,j] (let r be its best decision), D[i,l..r], and D[l..r,j].
When applying , due to the dependence from D[i,j] to D[i,j-1] and D[i+1,j], the δ-th frontier contains the states D[i,i+δ].
Hence, although it results in optimal parallelization to the standard sequential algorithm, the algorithm requires n-1 rounds and thus has O(nlog n) span. 
Achieving o(n) span may need new insights to redesign the dependencies. 














§ EXPERIMENTS


To demonstrate the practicability of our new algorithms,
we designed experiments for LCS and convex GLWS.
We implemented our parallel LCS algorithm and parallel GLWS algorithm in C++
using ParlayLib <cit.> to support fork-join parallelism and some parallel primitives (e.g., reduce).
Our tests use a 96-core (192 hyperthreads) machine with four Intel Xeon Gold 6252 CPUs and 1.5 TB of main memory.
Our code is available at <cit.>.

For parallel LCS, the existing parallel implementations we know of
<cit.>
cannot process inputs with 10^8 size,
so we compare our parallel LCS algorithm with the sequential version of our algorithm.
We test two random strings A[1..n] and B[1..n] with length n=10^8,
while controlling  (number of pairs (i,j) such that A[i]=B[j]) and k (the LCS length).
The pre-processing time to find all matching pairs is not counted into the running time.
<ref> shows the results when =10^8 and =10^9.
Our algorithm has up to 30× speed up than the sequential version.

For GLWS, we use the setting of the post office problem described in <ref>,
and compare our parallel algorithm with the sequential solution in <ref>.
We generate random data for n=10^8 and 10^9, and 

use different weight functions to control the output size k,
the number of post offices in the solution.
<ref> shows the result on different n and k.
The time for sequential algorithm does not change significantly,
because it has O(nlog n) work, independent of k.
For our algorithm, the running time varies with k due to the O(klog^2n) span. 
When k is small, our algorithm is 20× faster than the sequential algorithm
and achieves 30–40× self-relative speedup.





§ CONCLUSION AND FUTURE WORK


We systematically studied general approaches to parallelize classic sequential dynamic programming algorithms, 
particularly those with non-trivial optimizations such as decision monotonicity and sparsification.
We showed a novel framework, the , and apply it to different DP recurrences.
Theoretically, we gave the concept of optimal parallelism and perfect parallelism of a sequential algorithm, 
and showed that with a careful design, we can achieve optimal parallelism for the classic sequential DP algorithms in a (nearly) work-efficient manner, and perfect parallelism for some instances.
Practically, we show that our carefully-designed techniques do not include much overhead, and can outperform the original sequential version in a wide variety of cases.

We believe that the novelty of the techniques in this paper opens a list of interesting questions.
First, many of the new parallel algorithms are nearly work-efficient—we pick the most practical sequential algorithms
for each problem, but they can be off the best work bound by up to an O(log n) factor.
It is theoretically interesting to ask if we can match the best work bound in parallel. 
Second, among all these classic algorithms we looked at, one problem/algorithm that 
cannot be directly solved by the  is the RNA Secondary Structure <cit.>. 

Here, we may have 2^n different paths in a DP DAG, so applying  efficiently may need some complicated techniques.
Finally, we show how to faithfully parallelize the sequential DP algorithms. We are aware of other approaches <cit.> for LIS/LCS
that can achieve stronger worst-case span bounds using divide-and-conquer. 

Hence, an interesting direction is to see if we can redesign other DP algorithms in a similar form to achieve better worst-case span.








§ ACKNOWLEDGEMENT

This work is supported by NSF grants CCF-2103483, IIS-2227669, NSF CAREER Awards CCF-2238358 and CCF-2339310, the UCR Regents Faculty Development Award, and the Google Research Scholar Program.




ACM-Reference-Format









§ ADDITIONAL DETAILS FOR PARALLEL OPTIMAL ALPHABETIC TREES (OAT)


This section supplements additional details for the parallel OAT algorithm in <ref>.



 §.§ Sequential Algorithms for OAT


Hu-Tucker algorithm <cit.> was the first O(nlog n) algorithm for OAT,
and a simplified version Garsia-Wachs algorithm <cit.> was proposed later.
Both algorithms have two phases.
In the first phase, a certain tree called the l-tree (short for level-tree) is constructed from the input sequence.
Hu-Tucker algorithm and Garsia-Wachs algorithm differ in the manner to find the l-tree, but the final l-tree is the same.
In the second phase, the OAT can be constructed from the l-tree,
and each item in the OAT is at the same level as in the l-tree.

The key problem is how to compute the l-tree.
Here we describe the first phase of Garsia-Wachs,
as it is easier to understand than Hu-Tucker and the idea is used in the parallel algorithm by <cit.>.
Let the input weight sequence be a_1..n.
We denote a_i+a_i+1 be the i'th 2-sum, for 1≤ i < n.
A pair of consecutive elements (a_i,a_i+1) is said to be locally minimal if the i'th 2-sum is a local minimum in the sequence of 2-sums.

The Garsia-Wachs algorithm repeatedly performs the following steps until there is only one element in the sequence:

  
  * Find the left-most locally minimal pair (a_i,a_i+1).
  
  * Combine a_i and a_i+1.
  Make a new node x to be the parent of a_i and a_i+1 in the l-tree,
  and the weight of x is a_i+a_i+1.
  
  * Remove a_i and a_i+1 from the sequence.
  Insert x before the first element a_j where j>i and a_j≥ x.
  If such a_j does not exist, insert x at the end of the sequence.

Note that each newly generated node represents a tree, with both children as the two trees corresponding to the two nodes it merges in step (2).


At the end, there will be only one element in the sequence, which is the final output of the l-tree.

An l-tree can be easily converted to an OAT in parallel with O(n) work and O(log n) span, and we refer the audience to <cit.> for more details.



 §.§ The Parallel Algorithm by <cit.>


Here we focus on the first phase on how to construct the l-tree in parallel.

Larmore et al. observe that we can pick any locally minimal pair instead of the left-most one as in Garsia-Wachs,
which does not affect the resulting l-tree.
Hence, we can in parallel process all possible locally minimal pairs in one round.
However, in the worst case the number of rounds is still linear.
For example, if a_1..n is in increasing order, then only one locally minimal pair can be processed in one round.

To overcome this issue, Larmore et al. further proposed the concept of “valleys”.
A valley α=(a_i) is the largest contiguous subsegment of a that 1) contains a_i, and 2) contains no item larger than a_i.
Hence, a_i's valley contains the elements from a_i's subtree in the Cartesian tree of a.
Larmore et al. showed that several disjoint valleys in the sequence can be processed in parallel.
To understand this, let us consider a valley α.
Let Δ_α be the parent of valley α in the Cartesian tree.
Since α is maximal, the locally minimal pairs from this range will not interact with the outside elements.
Hence, we can consider valley α as an independent task, and repeatedly find and combine the locally minimal pairs inside α.
The only difference is that if the combined element x=a_i+a_i+1>Δ_α,
we mark it as x^* and put it in a separate queue.
After the parallel processing of disjoint valleys, we collect all marked elements and insert them into the right place.

In total, there can be at most n overlapping valleys.
Larmore et al.'s algorithm uses one special type of valley: the 1-valleys.
A 1-valley is defined to be a valley α such that for any node v in the subtree of α in the Cartesian tree,
if v has two children, at least one of the children is a leaf.
We can see that in the sequence the maximal 1-valleys are not overlapping, so we can process them in parallel.

They also proved that, in each round if we process all maximal 1-valleys and reinsert them as in the sequential order,
the number of maximal 1-valleys in the remaining sequence decreases by at least a half. 
Thus, the whole algorithm will finish in O(log n) rounds.

Before we show how to process a 1-valley, we introduce some definitions.
Let p(u) is the parent of u in the Cartesian tree of α.
We define a 1-valley α to be a regular valley if for u<v are two leaves in the Cartesian tree of α, then p(u)<p(v).
A sorted regular valley is defined as a regular valley with the minimum element in the first one.
Larmore et al. showed that a 1-valley can be transformed to a sorted regular valley in O(m) work and O(log m) span, where m is the size of this 1-valley.
A set of items S is defined to be lf-closed (short for leaf-father-closed) if for any a_i∈ S and a_i is a leaf in the Cartesian tree of α,
then p(a_i) is also in S.
The weight of S is defined as the sum of the weights of all items in it.
We define W_k as the minimum weighted lf-closed set with k items in it.
Larmore et al. showed that all W_k can be computed in O(mlog m) work and O(log m) span.

Now consider a sorted regular valley α with length m.
Our goal is to generate a sequence of forests _0..m'-1 (m'<m) so that each of them corresponds to a subset of subtrees in the l-tree.
Here 
_0 is empty.
The last forest _m-1, if m=m', corresponds to the l-tree if we merge all elements in this valley.
However, this process may not end here—as mentioned above, we will stop if the weights of subtrees exceed Δ_α.
Hence, we can have our final state _m'-1 with m'<m.
To compute _i, we will enumerate all _j for j<i, and find the best (minimum) transition from them.
Here a transition means to build an additional level in the l-tree, and the cost w(j,i)=W_2i-j.

Larmore et al. showed that the cost function w is convex, so computing _i is exactly a convex GLWS problem.
If we use <ref> to solve this convex GLWS problem, since each decision [i]=j will add another level to the trees in _i from _j, the effective depth is upper bounded by the overall l-tree height h.
Hence, the work and span for each 1-valley subproblem is O(mlog m) work and O(hlog^2 m) span, where m is the subproblem size.
Since the total size of 1-valleys in a recursive round is O(n), the work and span for one round are O(nlog n) work and O(hlog^2 n), respectively.
Multiplying this by O(log n), the number of recursive rounds, gives the cost bounds in <ref>.


Each _i (i>0) adds a new root p_i by combining two existing nodes (either a root from _i-1 or a single node from α).
Lawrence proves the following lemma:

  Let 0<i<m. Then there exists 0≤ j < i that:
  (1) the roots of _i are p_j+1..i;
  (2) the leaves of _i are W_2i-j;

Note that from _j to _i we added one level,
thus the weight of _j is the weight of _i plus the weight of W_2i-j.
To find the minimum weighted forests, we can reduce to a LWS problem on 0..m-1 by defining w(j,i)=W_2i-j.
Larmore et al. showed that the cost function w is convex, so computing _i is exactly a convex GLWS problem.
Since each decision [i]=j will add another level to the trees in _i from _j, 
After solving this convex LWS, the final forest can be recovered in O(mlog m) work and O(log m) span.
As <ref> shows, each [i]=j brings an extra level in _i.
So the height of the final forest is the same as the longest path in the dependency graph.
As mentioned above, the final forest is a subgraph in the final l-tree, and the OAT has the same height of the l-tree.
So we have k, the length of the longest dependency path, is upper bounded by the final OAT height h.



 §.§ Proof of Lemma <ref>


Here we provide the proof of Lemma <ref>. We first show the following lemma. 


  Let the weight of a subtree in an OAT as the total weight of all leaves in this subtree. 
  In an OAT, the subtree weight grows by at least twice for every three levels.



Here we denote p(v) as the parent node of v in the output OAT T and w(v) as the sum of the leaf weights in the subtree of node v.


We will show for any node a in the OAT, the great-grandfather of a (if exists) must have weight no less than 2w(a).
In the optimal alphabetic tree T, let node b be a's sibling and node c be a's parent.
WLOG we assume node a is the left child of node c.
Then, let d be c's sibling, e as c's parent, and f as e's parent.
Now the lemma is equivalent to w(f)≥ 2w(a).

We consider the first case when c is e's left child (case 1 in <ref>).
In this case, we must have w(d)≥ w(a), since otherwise a right rotation will decrement the total cost of the tree by w(d)-w(a), violating that T is an OAT.
Hence, w(f)≥ w(e)≥ w(d)+w(a)≥ 2w(a) in this case.

The second case is when c is e's right child.
There are two sub-cases.
First, if e is f's right child (case 2(a) in <ref>), then c, e, and f form another case 1.  Hence, we have w(f)≥ 2w(c)≥ 2w(a).
Second, if e is f's left child (see T_4 in case 2(b) in <ref>). We can double-rotate and get another valid alphabetic tree T_5.
As T_4 is the optimal alphabetic tree, we must have 2w(d)+3w(a)+3w(b)+w(g)≤ 2(w(d)+w(a)+w(b)+w(g)),
which leads to w(g)≥ w(a)+w(b).
So w(f)=w(d)+w(a)+w(b)+w(g)≥ w(a)+w(g)≥ 2w(a).

From all three cases above, we show w(f)≥ 2w(a), which proves the lemma that the weight doubles for every three steps up the OAT.


With integer weights in word size W, the weight of the root is at most O(W), and the weight of each leaf is at least 1.
In this case, the number of levels between them is at most O(log W). 
This proves <ref>.




Case 1: if c is e's left child (see T_1 in <ref>).
We can do a right rotate to get T_2, which is also a valid alphabetic tree.
As T_1 is the optimal alphabetic tree, we must have 2w(a)+2w(b)+w(d)≤ w(a)+2w(b)+2w(d), which leads to w(a)≤ w(d).
So w(e)=w(a)+w(b)+w(d)≥ 2w(a), and the parent of e (a's great-grandfather) must have weight at least 2w(a).

Case 2: if c is e's right child. There are two sub-cases.
(a) If e is f's right child, then c,e,f forms case 1 (symetric). So we have w(f)≥ 2w(c)≥ 2w(a).
(b) If e is f's left child (see T_4 in <ref>). We can do a double-rotate and get another valid alphabetic tree T_5.
As T_4 is the optimal alphabetic tree, we must have 2w(d)+3w(a)+3w(b)+w(g)≤ 2(w(d)+w(a)+w(b)+w(g)),
which leads to w(g)≥ w(a)+w(b).
So w(f)=w(d)+w(a)+w(b)+w(g)≥ w(a)+w(g)≥ 2w(a).

So the weight of the root of OAT is at least 2^h/3.
If W is the word size, we must have 2^h/3<W, which leads to h=O(log W).






