Fast Algorithms for $L_\infty$-constrained S-rectangular Robust MDPs
Authors: Bahram Behzadian
Marek Petrik
Chin Pang Ho
Venue: NeurIPS 2021 Poster (2021)
Keywords: reinforcement learning
robust Markov decision processing
rectangular ambiguity sets

Abstract: Robust Markov decision processes (RMDPs) are a useful building block of robust reinforcement learning algorithms but can be hard to solve. This paper proposes a fast, exact algorithm for computing the Bellman operator for S-rectangular robust Markov decision processes with $L_\infty$-constrained rectangular ambiguity sets. The algorithm combines a novel homotopy continuation method with a bisection method to solve S-rectangular ambiguity in quasi-linear time in the number of states and actions. The algorithm improves on the cubic time required by leading general linear programming methods. Our experimental results confirm the practical viability of our method and show that it outperforms a leading commercial optimization package by several orders of magnitude.

Fast Algorithms for Lo-Constraieed S-Rectangular Robust MDPs
Bahram Behzadian University of New Hampshire bahram@ cs. unh edu
Marek Petrik University of New Hampshire mpetrik@cs. unh edu
Chin Pang Ho
City University of Hong Kong clint. ho@ cityc. edu. hk
Abstract
Robust Markov decision processes (RMDPs) are a useful building block of robust reinforcement learning algorithms but be hard solve. This paper proposes a fast, exact algorithm for computing the Bellman operator for S-rectangular ro- bust Markov decision processes with -oonttrnrrd rectangular ambiguity sets. The algorithm combines a novel homotopy continuation method with bisection method to solve S-rectangular ambiguity in quasi-linear time in the number of states and actions. The algorithm improves on the cubic time required by leading general linear programming methods. Our experimental results confirm the prac- tical viability of our method and show that it outperforms leading commercial optimization package by several orders of magnitude.
Introduction
Markov decision processes (MDPs) are powerful framework for dynamic prciision-maiing prob- lems and reinforcement learning (Bertsekas and Tsitsiklis, 1996; Puterman, 2005; Sutton and Barto, 2018). The MDP model assumes that the exact transition probabilities and rewards available. However, these transition probabilities are typically unknown and must be estimated from sampled data. Such estimations are prone errors, and the MDP's solution sensitive the introduced statistical errors. In particular, the quality of the optimal policy degrades significantly even with small errors in the transition probabilities (Le Tallec, 2007).
Robust MDPs (RMDPs) mitigate MDPs' sensitivity to estimation errors by computing an optimal policy for the worst plausible realization of the transition probabilities. This set of plausible transi- tion probabilitics is known as the ambiguity set. this paper, we study RMDPs with S-rectangular ambiguity sets, which can be solved in polynomial time (Hansen, Miltersen, and Zwick, 2013). However, computing the worst-case realization of transition probabilities often requires solving a linear program (LP) another convex optimization problem. Modern solvers are powerful and efficient, but as the problem size grows, solving an LP for every state becomes computationally prohibitive (Ho, Petrik, and Wiesemann, 2018).
Most prior work has focused on RMDPs with L, constrriiee ambiguity sets because both conve- nient concentration inequalities (Weissman al., 2003; Petrik, Ghavamzadeh, and Chow, 2016; Russel, Gu, and Petrik, 2019) and fast algorithms (lyengar, 2005; Petrik and Subramanian, 2014; Ho, Petrik, and Wiesemann, 2020) exist for this scenario. The concentration inequalities play important role in the data-driven construction of high-confidence RMDPs. However, ambiguity sets defined by the L more natural and interpretable by human modelers (Givan, Leach, and
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

Dean, 2000; Delgado et al., 2016), and can significantly outperform L1 -bsed ambiguity sets in many cireumstances (Behzadian et al., 2021). Unfortunately, RMDPs with S S-rectangular ambiguity sets defined in terms of the Too ball can currently be solved only using which are complex and slow.
I.P solvers,
As our main contribution, we propose new, fast algorithm for solving RMDPs with L constrained ambiguity sets. Our algorithm combines new homotopy continuation method with a bisection method to achieve quasi-linear (SSA A log S) time complexity with respect to the number of states S and actions A. This computational complexity compares favorably with the cubic (((SA)) time complexity of general interior-poin LP algorithms. We identify new simplifying properties of the robust optimization problem defined over Lo balls develop our algorithms.
Although bisection and homotopy methods have been used previously in the context of robust MDPs, their use and assumptions differ significantly from this work. A bisection method was used to solve S SA-rectangular RMDPs (Nilim and EI Ghaoui, 2005), but their approach does not generalize to S-rectangular RMDPs that we target. Homotopy and bisection methods have been used to solve L Li-constrained ambiguity sets (Ho, Petrik, and Wiesemann, 2018, 2020) but these methods are based on sparsity properties of the Lı-norm, which do not hold for the -nom.. We elaborate on this important difference after we introduce our algorithm. The existing efficient algorithms are developed for the SA SA-rettangular RMDPs with L balls (Givan, Leach, and Dean, 2000), but they do not generalize to S-rectangular RMDPs. Developing fast optimization algorithms for S- rectangular RMDPs is especially challenging because optimal policies may need be randomized.
Several new, fast methods have been proposed recently for solving RMDPs more efficiently. They propose replacing the standard value and policy iteration methods with more efficient algorithms, such as forms of modified policy iteration (Kaufman and Schaefer, 2013; Ho, Petrik, and Wiesemann, 2020) or gradient descent (Grand-Clément and Kroer, 2021). Most of these accelerated methods can further benefit from the fast Bellman operator algorithms that we propose in this work.
The remainder of the paper is organized as follows: Section 2 describes the basic robust MDP frame- work. Then, Section 3 proposes a new homotopy method for solving SA-rectangular ambiguity sets, which serves building block for main contribution. In Section 4, we propose a bisection method that can solve, in combination with the homotopy method, RMDPs with S-rectangular am- biguity sets. Finally, Section 5 presents experimental results that show that our method is over 1000 times faster than using Gurobi, a leading commercial linear solver, when solving RMDPs with hundreds of states.
Notation: We reserve lower- and uppercase bold characters for vectors and matrices respectively. The symbol A denotes the probability simplex in IR. Finally, we use I, 1, 0 to denote an identity matrix, a vector of ones, and a vector of zeros respectively.
Preliminaries: Robust MDPs
This section surveys the basic properties of RMDPs; for example, please see (Iyengar, 2005; Wiese- mann, Kuhn, and Rustem, 2013; Ho, Petrik, and Wiesemann, 2020) for more details. We consider finite RMDP model with states S {1...., S} and actions A = {1......A A}. The agent takes an action a E A in state s E S, it receives a reward re.a E R. and it transitions to the next state s' E S with probability of Ps.a,st. The transition probabilities Ps.a,s' are unknown but restricted to be in an ambiguity set P C (ASSSX, The initial state is distributed according to Po A
We aim to compute a policy æ : + A from the set of stationary randomized policies II that max- imizes the expected d-dicouuntdd return p II x P + F for the worst-case transition probabilities:
max min p(T,P) well PEP
Here, p(z, p) is the standard discounted infinite-horizon return for a policy π defined as p(T, [E 3² r(S,, At) A, ~ x ((St),S4z ~ PPs..A, Sg ~ Po]. The optimization problem in (1) can be seen as a zero-sum game, where adversarial nature chooses transition probabilities from the ambiguity set order to minimize the agent's return. Since solving the general optimization problem in (1) is NP-hard (e.g., (Wiesemann, Kuhn, and Rustem, 2013)), most research has focused on RMDPs with S-rectangular and SA-rectangular ambiguity sets that can be solved in polynomial time (lyengar, 2005; Le Tallec, 2007; Wiesemann, Kuhn, and Rustem, 2013).

---

SA eA-rectangular ambiguity sets P are defined as Cartesian products of sets Ps.o S A for each state s and action a as P = {p E (A (A$)SxA Ps.a E Ps.a. 3 S, E A]. The intuitive interpretation of sf-rettangylrrttyy is that nature can choose the worst transition probabilities from sets Ps.a a for each state and action a independently. We focus on ambiguity sets bounded by -nom distance from nominal transition probabilities Pa.a E A defined as
Ps.a = {Pan E
||pr.a -P, S Kss0
(2)
where Ks.a 2 0 is the robustness budget, and the nominal transition probability Pa.a is typically estimated from samples of state transitions.
To streamline the definition of the robust Bellman operator, we follow the notation of Ho, Petrik, and Wiesemann (2018) and define a nature response function q IR+ R R that represents the nature's s response for a particular state 8, action a, and budget g as
(5,
min pEAs
{TH.. + ypupv Ps.a
P st)
(3)
Then, the SA-rectangular robust Bellman operator I
R$ R for a value function R is
(Tv) e max min
(E, a ESMA,e
The optimal value function v* E R must satisfy the robust Bellman optimality equation v* = To" and can be computed using value iteration, policy iteration, or other methods (lyengar, 2005; Ho, Petrik, and Wiesemann, 2020; Kaufman and Schaefer, 2013; Grand-Clément and Krver, 2021).
S-rectangular ambiguity sety relax the assumptions of SA-rectangular sets and compute less con- servative policies but with a higher computational complexity (Wiesemann, Kuhn, and Rustem, 2013). They are defined Cartesian products of sets Ps S (A$)4 for each state as P = {p (A (Paa)aEA E P.. Vs SJ. As with SA-rectangular sets, we also con- sider marginal ambiguity sets Pa defined in terms the -nor as
- (Pra)aEA E (AS)A
P. =
ΣE - Ps.alls S Kig aEAA
where 2 0 is the robustness budget and Ps.a is the nominal transition probability. The important distinction from the SA-rectangular setting is that Ks depends only on the state and not the action. The S-rectangular Bellman operator is then defined as
(Iu); = max min dEA* ESA
(5)
aeA
Notice that the S-rectangular Bellman operator allows for randomizing actions through the proba- bility distribution d, which improves robustness but introduces additional significant computational complexity (Wiesemann, Kuhn, and Rustem, 2013; Ho, Petrik, and Wiesemann, 2020).
The vast majority of RMDP methods employ value iteration and policy iteration principles and require computing the robust Bellman operator many times during their run (lyengar, 2005; Wiese- mann, Kuhn, and Rustem, 2013; Ho, Petrik, and Wiesemann, 2020). Therefore, is important that it can be computed efficiently. In the remainder of the paper, we develop new quasi-linear time algorithms for computing the robust Bellman operator.
3 Computing the SA-Rectangular Bellman Operator in Linear Time
In this section, we develop new quasi-linear time algorithm for computing the SA-rectangular robust Bellman operator defined by the Loc- norm. This entails solving the optimization in (4). The algorithm developed in this section also serves as the major building block of the S-rectangular algorithm described in Section 4. The remainder of the section is organized as follows: Section 3.1 first analyzes the I.P formulation of the function q, and, then, Section 3.2 uses these properties to develop a new, fast homotopy continuation algorithm.
3

---

1.00
Rolaxl O-funsiion: v
Index
Transtton probabiitty:
0.75
0.30
0
0.25
0.6
0.9
0.6 Size ofabbgguity
Size ambbluulty set
Figure 1: Function ((S) in Example 3.1.
Figure 2: Probabilities p (E) in Example 3.1.
Computing the SA-rectangular robust Bellman operator for fixed state s, action and value function requires one evaluate the nature response function (s,a(E, v) in (3). Because the sym- bols 8, a, are fixed throughout this section, we omit them in the notation. For example, we use g(E) instead of as.a v) and p place of Ps.a. To further eliminate clutter, let rssa" 11+19 +7- Then, the optimization problem in (3) can be formulated as the following parametric LP:
<(E)
min {p' z
p pll0 - =
PEA*
(6)
min pER
{
1'p = -
2 0, L....,,, s}.
The remainder this section develops fast algorithms for solving (6) for all values { 2 0.
3.1 Properties of Nature Response Function !
The LP in (6) can be solved using generic solvers, like Gurobi or Mosek, but these are impractically slow for solving RMDPs. The optimization in (6) can also be solved in quasi-linear time for any fixed { 2 0, as we summarize in Appendix C. The known quasi-linear algorithm is, unfortunately. insufficient for solving the S-rectangular robust Bellman operator in Section 4. In this section, we prove results that pave the way for solving (6) for all & 2 0 simultaneously in quasi-linear linear time, which enables efficient algorithms for both S- and SA-rectangular RMDPs.
will be convenient to use p" (E) to refer to an optimal solution in (6). To avoid unnecessary tech- nicalities, we assume that all elements of are distinct, which guarantees that the optimal solution (E) is unique. In practice, one may add an arbitrarily small value to the elements of z to ensure that they are all distinct.
To get some intuition into the form of the nature response function (()) and its optimal solution p" (E), consider the following simple example.
Example 3.1. Consider an RMDP with six states, action, (-1,0,,,,, T)T, and nominal transition probabilities
(0.0, 0.1, 0.3, 0.1, 0.3). The functions g(E) and (E) are depicted in Figures I and 2, where Figure 2 shows the evolution of each pı(E) using a different color for each i.
The following property of the function is indispensable for analysis and shows that g(E) is always of the form depicted in Figure 1. It follows from standard LP properties and is proved in Appendix A.I.
Lemma 3.2. The function q(E) is continuous, piecewise linear, non-increasing, and convex in &.
To develop an efficient algorithm, now analyze the structure of the bases of the LP (6). Recall that basis a subset linearly independent constraints in the LP, which must hold with equality. There are S constraints included in each basis because S is the number of optimization variables. Note that constraints may be active (or violated) without being included in the basis.
To represent basis in (6), we use sets Rg. NB. TB S (1....,S} S} to indicate which constraints are included the basis with their meanings summarized in Table If E Dn, we call it a donor, if E R.B, we call it receiver, and if E NNs we call it none. The set Tp = {1,..., S) Ra DB \NB represents the remaining indexes, and E TT called trader. Lemma 3.4 below justifies the names for these sets.

---

Index i e
Constraints in
RB (receiver)
Pi P: S E in
DB (donor)
P: p; S E in B
0.5
NB (none)
P: 2 0 in B
Figure 3: Composition of B for
ambiguit set
Figure 4: An illustration of Algorithm 1.
Our homotopy algorithm will leverage the specific behavior of the optimal solution (E) as func- tion of ç. Because each basis B represents a set of S linearly independent inequalities with S variables, there exists a unique solution (E) for any value Ç. Note that pn(() need not be optimal or feasible.
The following lemma establishes the properties of the bases in (6) that we need consider in our optimization; its proof can be found in Appendix A.1.
Lemma 3.3. Suppose that p* is optimal (6) for some & 2 Then, there exists a basis B such that (i) p" (E) (ii) sets RB. NB,NB, TB do not intersect, (iii) [R |Dal + [WB] TB S.
(iv) TB = 1, and (v) 2, < 2j < 2g, for each i E RB.j TB, DB UNB.
Lemma 3.3 is important because it limits the bases relevant to the optimization, which is crucial for building fast algorithms. In particular, shows that the sets R. D. partition the set S, and there is always exactly one trader. The lemma also shows that coefficients for receivers are smaller than the coefficient for the trader, which is smaller than the coefficients for donors and nones.
The following lemma establishes the rate of change of the linear function Pc((), which is the last necessary component for our homotopy algorithm. The lemma's proof is in Appendix A.I.
Lemma 3.4. The derivatives p VEPD(E) for any basis B that satisfies the properties in Lemma 3.3 are equal
PA = 1 if RB, p. p - DB. = = D= 1DBl =
RB if TB
for E S. Moreover, the slope is q = 4/de (B(() = ΣieRg ΣjED, Zj + E P7
Note that Lemma 3.3 shows that each i either a receiver, donor, a trader, or none. Lemma 3.4 then shows that with an increasing &. a donor donates its probability mass, receiver receives probability mass, a trader either donates or receives at variable rate, and a none remains unchanged.
3.2 Homotopy Algorithm
We are now ready to describe the proposed homotopy method and prove its correctness and complex- ity. Algorithm summarizes a conceptual version of the homotopy algorithm. As we discuss below, one needs to avoid computing the full gradient V Pp((() to achieve quasi-linear time complexity. The complete algorithm with quasi-linear runtime described in Algorithm 3 in Appendix A.2.
The main idea of Algorithm is simple: it iteratively computes the linear segments of q(E) for all where the optimal solution is Po = 5 with objective value go = p/z. Then, the algorithm tracks the optimal bases in ((E) as E increases. When the PB, (E) becomes infeasible with the increasing &, the algorithm finds continues until it arrives at basis with d/dg ((')) O; the function ş constant for all & > &'. Since ({{) is piecewise linear in (see Lemma 3.2), we obtain its full description from all optimal bases.
0. The algorithm starts with &
optimal basis B and
The following theorem proves the correctness of Algorithm 1. Informally, the theorem shows that the function q piecewise linear with breakpoints (points of non-linearity) only at & /. 1,..., Note that { = 1 because this is the upper bound on the Loomnorm of a difference of two discrete probability distributions, and, as result, the function q(E) constant for & > The proof can be found in Appendix A.I.
5

---

Algorithm 1 Homotopy method to compute g(E)
1: input: Objective and nominal probabilities p:
2: output: Breakpoints (E..... and (9:):=0.... such that q: =((E): :
3: Initialize Eo + 0, - 0. Po < p and do f ((5)) = pbz. TO = [S/2] and basis Bo such that:
4: Te. [ RB = {ili < Toh, DB. {ili Toh, No = {}
5: while di < 0 do
6: Compute maximum step size for B, to remain feasible (TB, (Ta)):
7: AC, + max (5 2
P. + PEP (Ex) 2 O, (P: + V&Pn, (5r) - p) E t7:
8:
Update breakpoints:
9:
P1+1 +- + A& V Pg, (Et); a++1 t p|+12: E-+1 < & + A
10:
Let Bt+1 next basis with the steepest slope (see Lemma 3.7 and Table 1);
11:
Let ++
end while
13:
Let ET+1 < and 4T + QT:
return: (Er)t=0,... +1, and (gt T+11
Type
B2+1
DB++1
R
TB
N
1: D + W
B1
DB, \ (2)
R.B.
TB.
NB, ( Nn, UTB. Nn.
2: T + N
B²
DB. DB, UTB.
RB. \ {m) RB, {m)
{m}
3: T D
B3
{m}
Table I: Possible types of basis change at breakpoint E++1 described in Lemma 3.7.
Theorem 3.5. Suppose that Algorithm 1 returns (E:) and (9))==0.....) T+1- Then g(a & + a1 El41) a ((E,) + a1- 1()++1) for [0, 1] and f O...., 11
We will refer to Figure 4 in order to provide the intuition that underlies the construction of Algo- rithm I and correctness. The figure depicts an example state of Algorithm at 2 and Line 10. The solid lines show the values QB, and 4B, when they are feasible and optimal. The dashed lines indicate when the bases are infeasible or suboptimal at each one of the breakpoints E2, The col- ored lines at & indicate the slopes for the possible candidates for B2. The algorithm chooses basis with the minimal slope.
The correctness of Algorithm 1 follows from the following three lemmas. The first lemma shows that the algorithm chooses the initial basis with the minimum possible slope.
Lemma 3.6. The basis Bo constructed in Line 3 of Algorithm steeper slope than any other basis that satisfies the conditions of Lemma 3.3:
feasible at E O and has
d/at qB, (0) S 4/ (0).
The second lemma shows that the next basis will be selected according to one of the rules in Table 1.
Lemma 3.7. Let basis B, be optimal for Ei+1 in Algorithm 1, such that p" (E1+1) = Ps, (E++1)
and q(Ex-1) qB. (E++1). Assume that PB. (E) is infeasible for & > E1+1. If B are all bases feasible for some { > ε,, then one with the steepest slope can be constructed as
B1 if
(P (E++1)); 0, for some E DD. (pv, (E(41)) = O, and TB. = [[T if PB, (E++1)) Ei+1, and TB. {7)
argmin BEB
((()))) 3
B
if
de
B3
where BB,BB² B3 are defined in Table I and E argmax Z
Lemma 3.7 shows that there are three possible types of basis change; any other possible choice of the basis would contradict the continuity of v(E) (Lemma 3.2). Recall also that Lemma 3.3 shows that there is always exactly one trader. The first type of basis change occurs when pt for a donor I E 'D reaches zero; the donor turns into none in the new basis. The second type of basis change

---

occurs when the trader probability mass becomes zero; the trader then turns into a none, and the receiver with the largest value becomes the new trader. The third type of basis change happens when the trader's gradient satisfies 4/dt P+( < -1 and its probability mass reaches to its lower bound for a given ç, making the basis infeasible for greater values of g. The trader then becomes a donor, and, again, the receiver with the largest value becomes the new trader.
Finally, the third lemma shows that the optimal basis Be identified at & remains feasible until Er+1- Note that the convexity of g() implies that the feasible basis also remains optimal.
Lemma 3.8. If B, is feasible and optimal at & in Algorithm 1, then it is also optimal on the interval [Er, & + AE.] computed in Line 6 of Algorithm 1.
We now turn to the computational complexity of Algorithm 3. As the following theorem shows, the number of iterations T in Algorithm is at most OSS). Unfortunately, keeping track of Pe in each iteration of Algorithm I requires also (SS) time leading to the overall time complexity of (S²)). To adapt Algorithm to run in quasi-lineds time, Algorithm 3, in Appendix A.2, generates the necessary values ε.: 9e without tracking the complete P. values. Its runtime is quasi-linear because it needs to sort the values of & to perform the optimization in Line 10 in constant time.
Theorem 3.9. Algorithm 1 terminates in at most O(S) iterations and can be adapted to run in OSS log S) time (see Algorithm 3 in Appendix A.2).
We conclude by discussing the relationship with the homotopy method proposed for solving RMDPs with the L1 ambiguity sets (Ho, Petrik, and Wiesemann, 2018). Although our algorithm also a homotopy method, it is based on analysis that departs significantly from earlier work. The simplify- ing properties for the Loe ambiguity sets differ significantly from the L -orm. When the ambiguity sets are defined as L1 balls, only two components of change at the time. Figure 2 illustrates that when the ambiguity sets are Loe balls, all components of p may change with the increasing Ç. The fast algorithm for the conntrriied RMDP relies on the more subtle structure of the optimal bases described in Lemma 3.3, which leads a more complex algorithm.
4 Computing the S-Rectangular Bellman Operator in Linear Time
In this section, we propose a fast algorithm for computing the robust Bellman operator (5) for S- rectangular RMDPs. We assume fixed state E S and omit the subscripts throughout the section. For instance, the nominal probabilities for state s and action a are denoted by P. E A4, We also assume fixed value function E R and let z, rs.a 1 + + v for A.
The fast algorithm for computing the S-rectangular robust Bellman operator builds on Algorithm 1. As Theorem 3.9 shows, the function ga defined in (3) is piecewise linear with ((S) linear segments that can be computed efficiently by Algorithm 3. Since qo is piecewise linear, it is easy to construct its inverse just by swapping &. and oi to get the following function:
aa (u) min {Ilp Pallx PEAS
p plz, S ku),
Va E A.
The function Un returns the budget that nature needs to achieve a response u. Using the function 9 we can reformulate the S-rectangular robust Bellman operator as:
Ed, . (a( (Iv), = max min dead dEAAgERA
Σ E. S
min 2
{
Σ Ia (u) S sn) m
(8)
WER
LAEA
aEA
aEA
The correctness of this formulation follows by standard duality arguments and is proved in Lemma A.3 in Appendix A.3.
The optimization in (8) is remarkable because its objective is me-dimensional function with one constraint. A natural algorithm to use with such an optimization problem is the bisection method outlined in Algorithm 2 (see Algorithm 4 in Appendix A.3 for more detailed algorithm). Algo- rithm 2 keeps an interval [Wminn Umax] such that the optimal u* satisfies that u² E [umnn Mmnax] In every time step, the algorithm bisects the interval [Umin: Mmax] in half and updates Umin, Umax in order to preserve that u² E [Mmin [m One may think of umin as the maximal known infeasible u in (8) and of Umax as the minimal known feasible in (8).
7

---

Algorithm 2 Bisection method for solving (7).
1: imput: Desired precision e, functions Ya Va A
2: output: à such that |u' - S 6, where u' is optimal in Equation (7)
3: Initialize bounds Minin A,
Umax max (Za)s:
4: while Umax Umin > 2ddo
5: Let + (Umin + Umax)/2;
6: if EMEA la (u) I then Umax + u else Umin < end if
7: end while
8: return: (Umin + 1/2
1,000
1,000
Relative Rurdin
Reletive Rurtime
100
200
10
100
300
400
100
200
300
400
Number of states
Number of states
Algorithm
Gurazl
Homorrayy
Algarihh
Gurbi
Bibection
Figure 5: Relative computation time (unitless) of our algorithms and an LP solver over nominal MDP in SA-rectangular (left) and S-rectangular (right) inventory management RMDPs.
The time complexity of Algorithm 2 depends the desired precision To remove this dependence on e, it is sufficient to replace the bisection by binary search over the breakpoints; we give the details of this method in Algorithm 4 in Appendix A.3. The following theorem, proved in Appendix A.3, summarizes the correctness and complexity of the proposed algorithms.
Theorem 4.1. The combined Algorithms and 2 compute the S-rectangular robust Bellman opera- for any state S and can be adapted (see Algorithms 3 and 4) run in time O(SA log S A)).
Numerical Results
This section compares the empirical runtime of Algorithms and with the runtime of Gurobi 9.1, leading LP solver. The results were generated on computer with an Intel i7-9700 CPU with 32 GB RAM; the algorithms are implemented in C++.
As the main benchmark problem, 2000). In this problem, the decision- maker must decide at every time step how much inventory to or- der. The number of states and actions in this problem corresponds to the holding capacity and order size respectively. This makes it easy to scale the number of states and actions and evaluate how the algorithms scale with problem size. To evaluate the performance of our methods we also consider the RiverSwim (RS) domain (Strehl and Littman, 2008) and the Machine Replace- ment (MR) domain (Delage and Ye, 2010). Please Appendix B for the detailed description of these domains.
the classic Inventory Management (IM) problem (Zipkin,
small problems,
Figure 5 shows the time to compute the robust Bellman operator for a single state in the inventory management domain. The x-axis represents the number states (maximum holding capacity) in the domain. The number of actions is the same the number of states. The y-axis represents the time to compute the robust Bellman operator divided by time to compute the standard (non-robust) Bellman operator. The results show that even in MDPs with few hundred states, the algorithms we propose are about 100 times faster than the leading LP solver. Interestingly, our algorithm an order of magnitude faster even for small problems. We use the robustness budget k 1.2, but the computation time is insensitive to the particular choice of k.
Table 2 compares the time to compute the robust policy for Machine Replacement (MR), RiverS wim (RS), and Inventory Management (IM) problems. The IM problem has 30 states. is worth empha- sizing that MR and RS are very small problems with less than 30 states, yet our algorithms are up to

---

Rect.
Algorithm
MR
RS
IM
Rect.
Algorithm
MR
RS
IM
SA
Algorithm Gurobi LP
<
3
SA
[Ho]-Alg.1
2
SA
2960
2240
9770
SA
Gurobi LP
92
363
1140
Algorithm 2
40
52
67
S
[Hg.-Alg.2
2
S
Gurobi LP
129
217
2740
S
Gurobi LP
79
317
2260
Table 2: Time (ms) to compute the robust pol- icy for S- and SA-rectangular RMDPs with Loo sets.
Table 3: Time (ms) compute the robust pol- icy for S- and SA-rectangular RMDPs with Ln sets (Ho, Petrik, and Wiesemann, 2020).
800 times faster than using an solver. This indicates not only that our methods scale well with the number of states but also that the constant overhead is quite small. For the sake of completeness, include in Table 3 the timing results obtained for the RMDP with L1 ambiguity sets. These results show that solving the connsttrre RMDP is more difficult than the In-constrained RMDP, but also that we can achieve similar dramatic speedups in Loe connttaiidd RMDPs as (Ho, Petrik, and Wiesemann, 2020).
Conclusion
We introduced new homotopy method for calculating robust Bellman operators for S- and SA- rectangular ambiguity sets constructed with -norm ball. Theoretically, we show that the worst- case time complexity of our algorithms is quasi-linear: O A log(S)). The algorithms also perform well in practice, outperforming a leading LP solver by several orders of magnitude.
In addition to being faster than sereraa-purpose LP solver, our algorithms are also much simpler. They make it possible to solve Loe connttraidd RMDPs without the cost and complexity of involving general LP solver. Although free and open-source LP solvers available, their performance falls significantly short of commercial ones. The algorithms we propose are also easy to combine with value function approximation methods in RMDPs (Tamar, Mannor, and Xu, 2014).
In terms of future work, we believe that is important to understand whether similar algorithms can be developed for RMDDs with more complex ambiguity sets, such defined using Wasserstein distance, Lormo or KL-divergence.
Acknowledgments
The authors would like thank Bence Cserna for discussions this topic and the reviewer for the comments that improved this paper. Partial support for the work was provided by the National Science Foundation (Grants IIS-1717368 and I1S-181,,75) the CityU Start-Up Grant (Project No. 9610481), the CityU Strategic Research Grant (Project No. 7005534), the National Natural Sci- ence Foundation of China (Project No. 72032005), and Chow Sang Sang Group Research Fund sponsored by Chow Sang Sang Holdings International Limited (Project No. 9229076). Any opin- ion, finding, conclusion, or recommendation expressed this material are those of the authors and do not necessarily reflect the views of the National Science Foundation and the National Natural Science Foundation of China.
References
Behzadian, B.; Russel, R. H.; Petrik, M.; and Ho, C. P. 2021. Optimizing Percentile Criterion using Robust MDPs. In International Conference on Artificial Intelligence and Statistics (AISTATS), 1009-1017.
Bertsekas, D. P.; and Tsitsiklis, J. N. 1996. Neuro- -ouro-dynamic programming. Athena Scientific, Bel- mont, MA.

---

Delage, E.; and Mannor, S. 2010. Percentile Optimization for Markov Decision Processes with Parameter Uncertainty. Operations Research, 58(1): 203-213.
Delage, E.; and Ye, Y. 2010. Distributionally robust optimization under moment uncertainty with application to data-driven problems. Operations Research, 58(3): 595-612.
Delgado, K. V.; De Barros, L. N.; Dias, D. B.; and Sanner, S. 2016. Real-time dynamic programming for Markov decision processes with imprecise probabilities. Artificial Intelligence, 230: 192-223.
Faísca, N. P.; Dua, V:; and Pistikopoulos, E. N. 2007. Multiparametric linear and quadratic pro- gramming, 1-23. John Wiley & Sons, Ltd.
Givan, R.; Leach, S.; and Dean, T. 2000. Bounded-parameter parameter Markov decision processes. Artificial Intelligence, 122(1): 71-109.
Grand-Clément J.; and Kroer, C. 2021. First-Order Methods for Wasserstein Distributionally Ro- bust MDP, arXiv:2009.0679))
Hansen, T. D.; Miltersen, P. B.: and Zwick, U. 2013, Strategy iteration is strongly polynomial for 2-player turn-based stochastic games with constant discount factor. Journal of the ACM (JACM), 60( I): 1-16.
Ho, C. P. Petrik, M. and Wiesemann, W. 2018. Fast Bellman Updates for Robust MDPs. In International Conference on Machine Learning (ICML), 1979-1988.
Ho, C. P. Petrik, M.; and Wiesemann, W. 2020. Partial Policy Iteration for L1-Robust Markov Decision Processes. arXiv:2006.09484.
Ibaraki, T:: and Katoh, N. 1988. Resource allocation problems: Algorithmic approaches. The MIT Press.
Iyengar, G. 2005. Robust dynamic programming. Mathematics of Operations Research, 30(2): 257-280.
Kaufman, L.; and Schaefer, A. 2013. Robust modified policy iteration. INFORMS Journal on Computing, 25(3): 396-410.
Le Tallec, Y. 2007. Robust, risk-sensitive, and data-driven control of Markov decision processes. Ph.D. thesis, MIT.
Nilim, A.; and El Ghaoui, L. 2005. Robust Control of Markov Decision Processes with Uncertain Transition Matrices. Operations Research, 53(5): 780-798.
Petrik, M.; Ghavamzadeh, M.; and Chow, Y. 2016. Safe policy improvement by minimizing robust baseline regret. Advances in Neural Information Processing Systems, 29: 2298-2306.
Petrik, M.; and Subramanian, D. 2014. RAAM The benefits robustness in approximating aggre-
gated MDPs in reinforcement learning. In Neural Information Processing Systems (NIPS).
Puterman, M. L. 2005. Markov decision processes: Discrete stochastic dynamic programming. John Wiley & Sons, Inc.
Rockafellar, R. T. 1996. Convex analysis. Princeton University Press.
Russel, R. H.; Gu, T.; and Petrik, M. 2019. Robust Exploration with Tight Bayesian Plausibility Sets. Multi- Dialii-Discppliaryy Conference on Reinforcement Learning and Decision Making (RLDM).
Strehl, A. L.; and Littman, M. 2008. An analysis of model-based interval estimation for Markov decision processes. Journal of Computer and System Sciences, 74(8): 3309-1331.
Sutton, R. S.; and Barto, A. G. 2018. Reinforcement learning: An introduction. MIT Press.
Tamar, A.; Mannor, S.; and Xu, 2014. Scaling Up Robust MDPs using Function Approximation. In International Conference of Machine Learning (ICML).
10

---

Weissman, T:; Ordentlich, E.; Seroussi, G.; Verdú, S.; and Weinberger, M. J. 2003, Inequalities for the LI Deviation of the Empirical Distribution.
Wiesemann, W.; Kuhn, D.; and Rustem, B. 2013. Robust Markov decision processes. Mathematics of Operations Research, 38(1): 153-183.
Zipkin, P. H. 2000. Foundations of inventory management. McGraw-Hil/rrwin ISBN 02561 13793.


Reviews:
[{"review_id": "U1LJHtwEHjb", "reviewer": "Reviewer_Xk7N", "summary": "The paper proposes a novel algorithm for computing the Bellman operator in robust MDPs with S-rectangular ambiguity sets on transition probabilities constrained by L-infinity norm. The algorithm is based on homotopy continuations and a bisection method, and results in orders of magnitude faster computation than existing algorithms based on linear programming, as verified empirically on problems with hundreds of states. This is a very impressive speed-up that allows significant increase in the size of robust MDPs that can be solved exactly. ", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The proposed solution based on homotopy continuations and bisection is very original and completely different from the existing method based on linear programming. Related work is adequately cited. \n\nThe quality of the paper is high. Theoretical analysis is sound, and the empirical evaluation is fair. The authors are comparing their method against the baseline linear programming method implemented using a fast commercial linear programming solver. The presented method is a complete piece of work. \n\nThe empirical demonstration demonstrates a spectacular speed-up with respect to the baseline method based on LP, but I am not so sure about the authors' claim that the proposed method has quasi-linear complexity in the size of the state space. The line graph for the computational time of the homotopy method in Fig. 5 looks very much like a straight line to me, which implies polynomial complexity. Furthermore, if I understand correctly, what is plotted in this graph is the ratio wrt the time it would take to do a Bellman back-up for a single state in regular, non-robust MDP, and if done directly, this operation is quadratic in the size of the state space for this inventory management example, where the number of actions is equal to the number of states. Wouldn't that make the absolute computational time of the homotopy algorithm scale not quasi-linearly with the size of the state space, but polynomially, with a rather high degree? I am not questioning the excellent computational performance of the algorithm, which alone looks sufficient for the acceptance of the paper to me, but perhaps the claim about its quasi-linear complexity should be explained a bit better.\n\nThe paper is very clearly written and organized. A minor typo, on line 139: \"one actions\" -> \"one action\". The results are very important, and of very high practical significance. The major computational speed-up provided by this algorithm allows solving of robust MDPs of sizes approaching industrial interest. The method does not appear to be very hard to implement, either. It is very likely that this method will be adopted by practitioners in the field of operations research.   \n\nUpdate: I am happy with the clarification the authors provided about the complexity of the algorithm, and trust it will find its way in the final version. I also realized that my rating was a bit of an outlier, so I revised to be more in line with the other ratings.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "I do not see any potential negative societal impacts of this work.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}, {"review_id": "KaWaPm2iL73", "reviewer": "Reviewer_7V4J", "summary": "The paper presents a fast algorithm for robust MDPs, specificially S-rectangular rmdps with L_inf constrained rectangular ambiguity sets. The main method proposed is to use homotopy and bisection. The authors compare their algorithm with Gurobi a leading LP solver on \"inventory management\" domain. ", "questions": "", "limitations": "", "rating": 5, "confidence": 1, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "First, let me tell you that I was not familiar with robust MDPs before reading this paper. \nThe paper is well written and it introduces the concepts used by referring to the appendix and / or other papers. Probably due to lack of my background knowledge, I had a very hard time following and understanding the paper. I carefully followed the definitions and equations at first, but understanding it fully required me to read other papers, and some chapters from the books. At some point, I gave up on completely understanding the concepts and read the paper without being able to verify the low-level details. In my opinion, a paper can be written in a way that makes it easy for people not familiar with the exact concept that the authors are building upon, but this is my personal opinion and I'm not an expert on robust MDPs, so I will leave this judgement to other reviewers.\n\nDespite my lack of knowledge, at the very high level, the paper appears to be sound and well written. Three very high level questions that I had were:\n - What is an example application of robust MDPs? The authors improve the speed upon LP solvers, but what kind of real world (or RL benchmark) problems we would be able to solve with this additional efficiency?\n - If the proposed algorithm is \"up to 800 times faster\" than LPs, would it be possible for the authors to apply it to problems with more than 30 states? Wouldn't that show the real strength of the algorithm?\n - Maybe I missed this information, but what stops this algorithm from working with RMDPs ambiguity sets such as L2 norm or KL divergeance?", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The paper does not comment on the potentials societal impact. It contains some information about technical limitations of the work as a future work (such as not handling RMDPs with more complex ambiguity sets).", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "CYJwNi3gumN", "reviewer": "Reviewer_Hc9Q", "summary": "This paper focuses on robust RL, especially developing a efficient algorithm for computing the Bellman operator for S-rectangular robust MDP (RMDP) with $L_{\\infty}$-constrained rectangular ambiguity sets. To achieve that goal, a homotopy continuation method combined with a bisection method is introduced. Besides, theoretical analysis is made to prove the correctness of the proposed method. Experimental results on inventory management problem, riverswim domain, and machine replacement domain are provided, comparing with \"Gurobi LP\" method.", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "- Originality: although this work extends the previous research on SA-rectangular RMDP, it is still novel to combine the homotopy method with bisection method, in order to solve RMDP with S-rectangular ambiguity sets. \n\n- Significance: robust RL is of great importance to RL community. This work focus on the aspect of computing an optimal policy for the worst plausible realization of the transition probabilities.\n\n- The relationships with the literature on computing the worst-case realization of transition probabilities are correctly reported; the theoretical analysis is detailed and sound; the empirical evaluation is demonstrated to support its claim on improving the efficiency.\n\n\n-------------------------------------------------------------------\nUpdates:\nI have read the author's response and other reviews. My concerns have been addressed by the rebuttal, therefore I am inclined to maintain my original score.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "- Notations and definitions should be clear. For exmaple, there is no definition for $\\xi$ in Eq.(3). Also, the subscript of $B$ in $\\mathcal{R},\\mathcal{D},\\mathcal{N},\\mathcal{T}$ in line 149 should denote the basis, but there is no definition.\n\n- \"Fig\" and \"Figure\" are both used, but it should be consistent to use either one only.\n\n- In line 135, it is reported that \"an arbitrarily small value\" will be added to the elements of $z$ to \"ensure that they are all distinct\". It's better to elaborate this a little bit in the appendix.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "7eTpt0vKPE", "reviewer": "Reviewer_7bXv", "summary": "This paper proposes to solve robust Markov decision processes, which are building blocks of robust reinforcement learning. The article derives new properties of S- and SA-rectangular Bellman operators. Based on these properties, the proposed algorithms for S- and SA-rectangular Bellman operators are exact and efficient in practice. Moreover, they compare favorably with the state-of-the-art, the best linear program solvers on the market. ", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "** Originality \n\nTo the best of my knowledge, the contributions of the paper are original. The authors identify new structural properties in the SA-rectangular Bellman operator, show how the homology algorithm can exploit them, and then extend the approach to the case of the S-rectangular  Bellman operator. That is quite a remarkable technical work. Furthermore, the complexity analysis of the algorithm proves a significant reduction compared to state-of-the-art, and the empirical evaluation confirms it. \n\n** Quality \n\nI did not find any false notes in the derivations. However, I did not do a detailed reading of the evidence. The empirical analysis was restricted to the performance evaluation of the Bellman operators. That is already good. For a better assessment of the contribution, performances as soon as the operators are integrated into the algorithms for solving robust Markov processes would have been welcome.  \n\n** Clarity \n\nThe paper is very well written and structured. Furthermore, the preliminaries of the paper allow the reader not familiar with the subject to have a good understanding of the problem attacked in this paper. \n\nExample 3.1 deserves a little more elaboration. I failed to interpret Figure 2 correctly. \n\nLine 139 - \"one actions\" -> one action? \n\n\n\n** Significance \n\nThe work is essential to the robust reinforcement learning community and may also interest the broader NeurIPS community. The impact is on the efficiency of solving the Bellman S- and SA-rectangular operator. \n\nThe significance would have been even greater if the authors had demonstrated the impact of the developed operators into algorithms for solving robust Markov decision processes. \n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "N/A", "review_text": "", "extra_scores": {"ethical_concerns": "N/A", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}]

Author rebuttal:

Meta-review: This appears to be a solid paper with the most confident reviewers having the most favorable opinions. There reviewers had some minor quibbles that should be easily addressed in the final version.

Decision: Accept (Poster)