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 -oonttrntrd 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 Fe.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-funstion: 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)
MER
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 pereraa-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": "n949zB0VSd", "reviewer": "Reviewer_hdeg", "summary": "In this work, authors consider S-rectangular roust MDP with L infinity constrained rectangular ambiguity and suggest a fast algorithm for computing the Bellman operator for such MDPs. Ambiguity sets are considered in the paper defined through the L-infinity norm. Specifically, in the definition, the distance between transition probability p_{s,a} and the nominal transition probability p_{s,a} is defined through L_{infinity} norm. In the prior work (Ho, Petrik and Wiesemann, 2018) instead of L_{\\infinity} norm authors used L_1 norm and similarly proposed a homotopy algorithm.\n", "questions": "", "limitations": "", "rating": 5, "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": "I found this paper a really interesting read however, I have several concerns concerning the novelty. It is not clear to me what are the major differences to the work in (Ho, Petrik and Wiesemann, 2018). Can the authors please clarify the main contribution beyond replacing the L_1 norm from (Ho, Petrik and Wiesemann, 2018) with a L_infity norm? Does the analysis follow from (Ho, Petrik and Wiesemann, 2018) or new innovations were needed to propose the current algorithm? If so, can the authors please elaborate on what are those main differences? \n\n1. It has not been clear to me what is the benefit (in terms of qualities of the solution when to using L_{infinity} norm instead of L_1 norm in the definition of ambiguity set? I think there has been an attempt at motivation but it was still not clear to me. Can the authors in their rebuttal clarify this through an example maybe? \n\n2. We know that  L_{infinity} and L_1 norms are equivalent (up to a constant). Couldn't we have applied the homotopy algorithm from (Ho, Petrik and Wiesemann, 2018) to ambiguity sets defined via L_{infinity} (by appropriately scaling robustness budget \\kappa_s)? If not, could the authors please clarify why not?\n\n3) The distinction between prior work (Ho, Petrik and Wiesemann, 2018) provided in lines 228-235 was not very convincing. Especially from the point of practicality of quality of the final solution. Can the authors please help me understand this point as well?\n\n4) On the numerical side, can we compare the authors please compared to the solution versus method in (Ho, Petrik and Wiesemann, 2018) directly in one table? Reading that part of the paper I was confused about what it is meant to say? Does this paper recover the old results, outperforms them, or is the comparison all together isn't valid to conduct? For instance, Table 2,3 does not show their relative performance, but only shows how they stand against Gurobi solver.\n\nMinor Typos: \ntypo: line 24 - Algorithm 1, not Algorithm 3.\ndefinition of \\boldsymbol{z}_a uses r_{s,a}, but previously you decided to omit subscripts s.\n\n\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Limitations were not clear to me from reading the paper especially compared to previous SOTA.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "TiAXc7cXw9", "reviewer": "Reviewer_So7E", "summary": "This paper proposes a new algorithm that can provide fast and exact solutions for the Bellman operator for S-rectangular robust Markov decision processes with L_\\infty-constrained rectangular ambiguity sets. The authors show that the solving time gain can be a cubic factor in comparing to leading general linear programming methods. This insight is further demonstrated through experimental results in comparisons to those general-purpose LP solvers. \n", "questions": "", "limitations": "", "rating": 6, "confidence": 2, "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": "In overall, I find the paper delivered very clear motivations on why solving RMDP with L_\\infty-constrained rectangular ambiguity sets is needed and challenging. It's very convincing that the proposed algorithm has come with theoretical results.\n\nThe experiment results show a lot of promising. The proposed algorithm seems to outperform the general-purpose LP solvers in terms of run time. However, given this runtime advantage, it might be also more interesting and convincing if there are results on problems with larger state spaces, at least similar to the tasks used in (Ho, Petrik, and Wiesemann, 2020) (discretized Cart-pole). \n\n\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "KWwKSN_6r8", "reviewer": "Reviewer_n2xU", "summary": "The authors in the paper present a new, fast algorithm for solving RMDPs with L∞-constrained ambiguity sets. In particular, they propose a new homotopy method for solving SA-rectangular ambiguity sets and a bisection method that can solve, in combination with the homotopy method, RMDPs with S-rectangular ambiguity sets. The experimental results show that our method is over 1, 000 times faster than using Gurobi, a leading commercial linear solver when solving RMDPs with hundreds of states. Overall, this paper proposes interesting methods to faster Bellman operator algorithms, which leads to solving RMDPs more efficiently.\n \n", "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": " Major:\n*All theorems and lemmas are well-established for the new, fast algorithm for solving RMDPs with L∞-constrained ambiguity sets. \n*The algorithms are well-developed and well-described. \n*The structure of the method is clear-represented, and the paper is well-structured. Lemma 3.3 limits the bases relevant to the optimization, which is crucial for building fast algorithms. \n*The fast algorithm for the L∞-constrained RMDP relies on the more subtle structure of the optimal bases described in Lemma 3.3, which leads to a more complex algorithm. \n*The experiment proves that the algorithm is superior (simpler and faster) to a general-purpose LP solver.\n*One drawback is that the paper does not discuss the non-rectangular complex ambiguity set case.\n \nMinor:\n \n*Line 58: …the basic Robust MDP framework.\n*X label in figure 4: the size of ambiguity set -> size of ambiguity set (missing the tail of letter g while cutting the image).\n*Line 225: optimizaton -> optimization.\n*Line 237: a fast algorithm for compute… -> a fast algorithm for computing….\n \nLimitations and societal impact:\nThe authors introduced a new, fast method for calculating robust Bellman operator for S- and SA- rectangular ambiguity sets constructed with L∞-norm ball. The algorithms are easy to combine with value function approximation methods in RMDPs. Whether similar algorithms can be developed for RMDPs with more complex ambiguity sets is a valid direction for future work.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "N/A.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "10 hrs"}, "extra_text": {}}, {"review_id": "EgO-MCPd5Ua", "reviewer": "Reviewer_cHzA", "summary": "This work studies the Robust MDPs (RMDPs) with S-rectangular ambiguity sets. The authors propose a fast algorithm to compute the robust Bellman operator over S-rectangular $L_\\infty$-ball ambiguity sets, which achieves quasi-linear $\\mathcal{O}(SA \\log S)$ time complexity.", "questions": "", "limitations": "", "rating": 6, "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: \nIn contrast with the RMDPs over $L_1$-ball ambiguity sets, the authors study the RMDPs with $L_\\infty$-constrained ambiguity sets. Therefore, they propose two algorithms to solve S- and SA-rectangular RMDPs over $L_\\infty$-ball ambiguity sets, which differ from the existing methods for solving RMDPs with $L_1$-constrained ambiguity sets.\n\nQuality: \nIn Lines 280-281, the authors claim that \"the computation time is insensitive to the particular choice of $\\kappa$\". However, they did not provide evidence to support the claim.\n\nClarity: \nThis paper is well written.\n\nSignificance: \nThis work studies the RMDPs over S-rectangular $L_\\infty$-ball ambiguity sets, which are more natural and interpretable than those over $L_1$-ball ambiguity sets. Experiments demonstrate that the proposed method outperforms the leading LP solver in terms of efficiency.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "15 hours."}, "extra_text": {}}]

Author rebuttal:

Meta-review: The reviewers and AC discussed the paper. The author response was helpful in clarifying some confusion. There was agreement that the paper is interesting. There are still some concerns about novelty and missing technical details, so this should be further clarified in the final version (in particular, the differences with Ho, Petrik and Wiesemann, 2018). 
Decision: Accept (Poster)