Ordinal Optimization Through Multi-objective Reformulation

We analyze combinatorial optimization problems with ordinal, i.e., non-additive, objective functions that assign categories (like good, medium and bad) rather than cost coefficients to the elements of feasible solutions. We review different optimality concepts for ordinal optimization problems and discuss their similarities and differences. We then focus on two prevalent optimality concepts that are shown to be equivalent. Our main result is a bijective linear transformation that transforms ordinal optimization problems to associated standard multi-objective optimization problems with binary cost coefficients. Since this transformation preserves all properties of the underlying problem, problem-specific solution methods remain applicable. A prominent example is dynamic programming and Bellman's principle of optimality, that can be applied, e.g., to ordinal shortest path and ordinal knapsack problems. We extend our results to multi-objective optimization problems that combine ordinal and real-valued objective functions.

Introduction
Ordinal objective functions occur whenever it is only possible to represent the quality of an element by a ordered category and not by a numerical value. As an example, consider the problem of finding optimized routes for cyclists in a road network: While edges may be associated with different categories like asphalt, gravel or sand-or, when related to safety considerations, very safe (there is a bicycle path), neutral (a quiet road) or unsafe (a main road without bicycle path)-such categories do not immediately translate into monetary or cost values. Bi-objective shortest path problems with route safety criteria are addressed, for example, in the web application geovelo and in the associated publications Kergosien et al. (2021); Sauvanet & Néron (2010). In these references, only two categories are considered (safe or unsafe edges), and the safety criterion is translated into a cost function that evaluates the total length of unsafe route segments. In contrast, an ordinal shortest path problem is investigated in Schäfer et al. (2020). A major difficulty when considering ordinal objective functions is that "optimality" may be defined in many different ways. Schäfer et al. (2020) suggests an optimality concept that is based on sorted category vectors. A similar concept is used in Klamroth et al. (2021), where matroid optimization problems with one real-valued and one ordinal objective function are investigated. A different perspective is proposed in Schäfer et al. (2021) who define ordinal optimality for knapsack problems on the basis of numerical representations for the categories.
In this paper, we consider general combinatorial optimization problems and provide a new cone-based interpretation of the optimiality concept for ordinal objectives suggested in Schäfer et al. (2021). In particular, we interrelate ordinal optimality with the classical concept of Pareto optimality for an associated multi-objective optimization problem. For a general introduction to multi-objective optimization we refer to Ehrgott (2005). Since the underlying transformation of the objective function is linear and bijective and hence preserves the combinatorial structure of the respective problems, our results immediately lead to efficient solution strategies as, for example, dynamic programming for shortest path problems. The respective transformations are based on a representation of dominance relations by cones.
The paper is organized as follows. In Section 2 we review three optimality concepts for ordinal optimization based on Schäfer et al. (2021). These results and optimality concepts are then extended and re-interpreted as a special case of cone-optimality (see, e.g., Engau, 2007) in Section 3. A detailed analysis of the properties of the corresponding ordering cones then leads to a linear transformation of ordinal optimization problems to associated multi-objective optimization problems with the Pareto cone defining dominance and optimality. This transformation is used in Section 4 to formulate a general algorithm for solving ordinal optimization problems. Furthermore, we investigate the relation between the definition of ordinal optimality and the weight space decomposition of the associated multi-objective problem. In Section 5 we extend our results to more general problem types with additional real-valued objective functions. We conclude in Section 6 with a summary and an outlook on future research.

Problem Definition
We consider combinatorial optimization problems with an ordinal objective function. In general, an ordinal optimization problem (OOP) can be formulated as " ordinally minimize " o(x) where X is the set of feasible solutions. We assume that X is a subset of the power set of a discrete set S, i.e., X ⊆ 2 S . Every element of S is assigned to one of K ordered categories. This assignment is encoded by a mapping o : S → C with C = {η 1 , . . . , η K }. We assume that category η i with i ∈ {1, . . . , K − 1} is strictly preferred over category η i+1 , written as η i ≺ η i+1 . The objective function of a feasible solution x = {e 1 , . . . , e n } is given by the ordinal vector o(x) = sort(o(e 1 ), . . . , o(e n )), where the operator sort( ) means that the components of o(x) are sorted w.r.t. non-decreasing preferences, i.e., o 1 (x) ≺ = o 2 (x) ≺ = · · · ≺ = o n (x). Note that different feasible solutions may have different numbers of elements, and hence the length of the ordinal vector o(x) may vary for different x ∈ X.
Instead of using the un-aggregated, ordered ordinal vector o(x) one can count the number of elements in a feasible solution per category. Accordingly, we use the counting vector c : X → Z K ≧ with Z K ≧ := {y ∈ Z K : y i ≥ 0 for all i = 1, . . . , K}. Thereby, the i-th component of c(x) equals the number of elements in x which are in category η i , i. e., c i (x) = |{e ∈ x : o(e) = η i }|. Obviously, there is a one to one correspondence between the vectors o ∈ R n and c ∈ R K , since the ordinal vector o can be determined from a given counting vector c by o i (x) = η j with j = argmin{j ∈ {1, . . . , K} : i ≤ j l=1 c l (x)}. Note again that the number of elements n of a feasible solution, and hence the length of the ordinal vectors o ∈ R n , may vary while the number of categories K and therefore the length of the counting vectors c ∈ R K is fixed. Hence, we get the following formulation of an ordinal counting optimization problem (OCOP) which will be shown to be equivalent to problem (OOP) for an appropriate definition of "ordinal minimization". In the following, we also consider an incremental tail counting vectorc ∈ R K that counts, in its i-th component, the number of elements of a feasible solution x which are in category η i or worse, i.e.,c i (x) = |{e ∈ x : η i ≺ = o(e)}| = K j=i c j (x). In particular, the total number of elements of a solution x is given in the first component ofc, i.e., |x| =c 1 (x) = K i=1 c i (x). As an example, consider the shortest path problem shown in Figure 1 together with the outcome vectors o(x) (for problem (OOP)) and c(x) (for problem (OCOP)) for all feasible solutions x ∈ X. In addition, the incremental tail counting vectorc(x) is given for all x ∈ X.

Optimality Concepts for Ordinal Objective Functions
In the following, we review three different concepts of optimality for ordinal optimization. All of them try to answer the question what minimization could mean for the problems (OOP) and (OCOP). The first concept is to use numerical representations that assign a numerical value to every category such that the order of the categories is respected. In this context, a numerical representation respects the order of the categories whenever the numerical value of a o cc e 4 e 5 e 6 e 7 e 8 Figure 1: Instance of an ordinal shortest path problem. A dotted-green edge is in the best category η 1 , a dashed-orange edge is in category η 2 and a solid-red edge is in the worst category η 3 . All possible s − t paths x i , i = 1, . . . , 6 and their respective objective function vectors are given.
better category is strictly smaller than the numerical value of a less preferred category. If we take the sum over all numerical values of a vector o(x ′ ) for a feasible solution x ′ , we obtain a unique numerical value that can be compared to the corresponding numerical value of another feasible solutionx. A feasible solution x ′ is called efficient if there is no other feasible solutionx which is better w.r.t. all numerical representations. The second concept is to maximize the number of elements in the good categories, and the third concept is to minimize the number of elements in the bad categories. After the formal introduction of these three optimality concepts, we investigate their interrelation.
Optimality by Numerical Representations. We consider the problems (OOP) and (OCOP). The concept of optimality by numerical representation for ordinal objectives as introduced in Schäfer et al. (2021) is based on a previous and more general work of Fishburn (1999). It assigns an order preserving numerical value to each category. Following Schäfer et al. (2021), we call a function ν : Note that we assume strictly ordered categories, i.e., there are no categories that are indifferent. As a consequence, we do not allow ν(η i ) = ν(η j ) for i = j since this would make two different categories indistinguishable in the numerical representation. Let V denote the set of all numerical representations for a given number of categories K.
For a given numerical representation ν, we define the numerical value of a feasible solution x = {e 1 , . . . , e n } ∈ X w.r.t. ν (cf. Schäfer et al., 2021) as The numerical value ν(x) of a feasible solution x ∈ X can be evaluated in different ways by re-arranging the terms and using the counting vector c or the incremental tail counting vectorc, respectively: An illustration for different ways to evaluate ν(x) is given in Figure 2.
Example 1. We apply the concept of numerical representations to the shortest path problem given in Figure 1. As a motivation, suppose that there are two decision makers A and B who have to select a most preferred path. They would agree, for example, that path x 1 is worse than path x 2 , because x 1 has a dashedorange edge more than x 2 and, other than that, their outcome vectors are the same. But they do not agree on the question whether x 2 or x 5 is preferred, because decision maker A chooses the numerical representation ν A (η 1 ) = 1, ν A (η 2 ) = 2 and ν A (η 3 ) = 5, while decision maker B chooses ν B (η 1 ) = 2, ν B (η 2 ) = 3 and ν B (η 3 ) = 4. Therefore, decision maker A would prefer path x 5 because ν A (x 5 ) = 4 < 6 = ν A (x 2 ) while decision maker B would prefer path x 2 because ν B (x 2 ) = 6 < 7 = ν B (x 5 ). Hence, the path x 2 does not ordinally dominate the path x 5 , i.e., x 2 is not better than x 5 for all numerical representations. Similarly, the path x 5 does not ordinally dominate the path x 2 since also x 5 is not better than x 2 for all numerical representations. Figure 2: Consider an example with n = 7 and K = 4. Different ways to compute the numerical value of a feasible solution x with o(x) = (η 1 , η 2 , η 2 , η 3 , η 4 , η 4 , η 4 ) ⊤ , c(x) = (1, 2, 1, 3) ⊤ , ν(η 1 ) = 2, ν(η 2 ) = 4, ν(η 3 ) = 7 and ν(η 4 ) = 8 are illustrated. The different colors represent the summands and visualize the different slicing strategies.
Definition 2 (cf. Schäfer et al., 2021). Let x ′ ,x ∈ X be feasible solutions. Then, and c(x ′ ) weakly ordinally dominates c(x), denoted by , respectively, if and only if for every ν ∈ V, it holds that , respectively, if and only if x ′ weakly ordinally dominatesx and there exists ν * ∈ V such that ν * (x ′ ) < ν * (x). 3. x * ∈ X is called ordinally efficient, if there does not exist an x ∈ X such that x x * . 4. o(x * ) and c(x * ) are called ordinally non-dominated outcome vectors of Problem (OOP) and (OCOP), respectively, if x * is ordinally efficient.
Optimality by Maximization of Elements in Good Categories. Another optimality concept in ordinal optimization is to maximize the number of elements in good categories. The intuition behind this concept is that solutions with many good elements are to be preferred over solutions with few good elements. The drawback, however, is that this concept rewards solutions with larger numbers of elements as long as these are in (relatively) good categories, which may not be wanted in practice. This optimality concept is defined only for the problem (OCOP), as we need the counting vector c for its definition.
Definition 3. We say x ′ weakly head-dominatesx, denoted by x ′ ≧ hx or Optimality by Minimization of Elements in Bad Categories. The drawback that longer solutions may be preferred over shorter solutions, as long as the elements are in good categories, can be avoided by taking the converse perspective, i.e., when minimizing the number of elements in the bad categories. Again, this optimality concept is defined only for the problem (OCOP).
Definition 4. We say x ′ weakly tail-dominatesx, denoted by (2) , if and only if (2) holds and c(x ′ ) = c(x). Moreover, x * ∈ X is called tail-efficient if there is no x ∈ X such that x t x * . The corresponding outcome vector c(x * ) is called tail-non-dominated.
Remark 5. Note that head-dominance as well as tail-dominance are equivalently defined on the feasible set X ⊆ 2 S and on its image set c(X) ⊆ R K . The definitions immediately extends to the complete R K .

Properties of and Interrelations between Optimality Concepts for Ordinal Optimization
In addition to the concepts described above, there are further ways to define efficiency. This has been done, for example, in Schäfer et al. (2020) for the ordinal shortest path problem. Their definition has the disadvantage that Bellman's principle of optimality (see Bellman, 1957) does not hold in general, i.e., not every subpath of an efficient path is necessarily efficient w.r.t. this optimality concept. The definition of head-optimality has the same disadvantage, see Remark 10 below for more details. In contrast, the definitions of ordinal optimality and tail-optimality can be proven to be equivalent. Moreover, they are compliant with Bellman's principle of optimality. Note that, for the special case of a knapsack problem, this was shown in Schäfer et al. (2021).
Proof. First let ν ∈ V be an arbitrary numerical representation. Then ν( Note that since n < m, and due to the above assumptions, there must exist a category η j , j ∈ {1, . . . , K} such that c j (x) < c j (x ′ ). Therefore, there exists a numerical representation ν * such that ν * (η j ) > 0 and thus ν * (x) < ν * (x ′ ), which concludes the proof.
Note that the condition of Lemma 6 is always satisfied ifx, x ′ ∈ X withx x ′ . Thus, in Example 1 the path x 2 is always preferred over the path x 1 . This is also the case for tail-efficiency, but not for head-efficiency, see Remark 10 below. In many application contexts this is meaningful property, since adding additional elements to a solution (no matter from which category) does generally not improve the solution quality. As an example, we refer again to paths representing bicycle routes as in the app geovelo, where we are not interested in routes that are unnecessarily long.
Lemma 7 (Schäfer et al. 2021). The ordinal dominance relation ≺ = defined on the feasible set X is a preorder, i.e., it is reflexive and transitive.
for every ν ∈ V by definition and therefore, we have shown transitivity.
Note that the ordinal dominance relation ≺ = is in general not antisymmetric on the feasible set X since two different feasible solutions may have the same number of elements in each category like the paths x 3 and x 4 in Example 1.
The following two results show that (weak) ordinal dominance and (weak) tail-dominance are actually equivalent on the feasible set X.
Proof. The proof is a simplified variant of the proof in Schäfer et al. (2021). Note, that they consider maximization problems while we consider minimization problems.
First we show by contradiction that The idea of the proof is to make the bad categories η j * , . . . , η K very expensive, such that an element of this category can not be replaced by elements of the lower categories. Hence, we define the numerical representation For the other direction we use the reformulation of ν(x), which is visualized in Figure 2(c). It follows that for any ν ∈ V

The inequality holds because of the assumptionc
For the other direction it is sufficient to show that c( Now the result follows analogously to the proof of Lemma 8 with exchanged roles of x ′ andx. Remark 10. Lemma 8 (and thus also Lemma 9) does not hold in general for the relation h . As a counter example, consider the paths x 1 and x 2 with c(x 1 ) = (1, 1, 1) ⊤ and c(x 2 ) = (1, 0, 1) ⊤ from Figure 1. Obviously, x 1 headdominates x 2 . But for every numerical representation it follows ν(x 1 ) > ν(x 2 ), which contradicts Note that the crucial point in the counter example given in Remark 10 is the different cardinality of the solutions, |x 1 | = |x 2 |. For ordinal optimization problems with fixed cardinality, for which matroid optimization problems as studied in Klamroth et al. (2021) are an example, it can be shown that head-and taildominance are equivalent.
Lemma 11. If all feasible solutions have the same cardinality, i.e., if |x ′ | = |x| for all x ′ ,x ∈ X, then head-and tail-dominance as defined in Definitions 3 and 4, respectively, are equivalent.
Lemma 12. The relation ≦ t is a partial order on R K , i.e., it is reflexive, transitive and antisymmetric. Moreover, the relation t is a strict partial order on R K , i.e., it is irreflexive and transitive.
This implies that u = v and hence ≦ t is antisymmetric. Therefore, ≦ t is a partial order. Now consider the relation t . Since u t u for all u ∈ R K , it holds that t is irreflexive. It remains to show that t is transitive. Towards this end, consider three vectors u, v, w ∈ R K such that u t v and v t w. This implies . . , K and u = w, i.e., u t w. Consequently, we have shown that t is irreflexive and transitive which concludes the proof.
As a consequence of the discussion in Sections 2.2 and 2.3, we focus in the following on the ordinal optimization problem (OOP) w.r.t. optimality by numerical representation, or equivalently, on the ordinal counting optimization problem (OCOP) w.r.t. the concept of tail-dominance.

Ordinal Optimality versus Pareto Optimality: An Interpretation Based on Ordering Cones
A prominent example of an order relation in the context of multi-objective optimization is the component-wise order or Pareto order, see, e.g., Ehrgott (2005). For two vectors u, v ∈ R K , we write and say that u Pareto dominates v if and only if u v. The associated weak and strict component-wise orders, respectively, are defined by Note that ≦ defines a partial order in R K while as well as < define strict partial orders in R K .
In the following, basic properties of orders and their relation to cones are discussed and illustrated at the example of the Pareto order. This leads to a novel perspective on ordinal optimality in comparison to Pareto optimality.

Orders and Cones
Orders and cones are closely related. The following review of basic concepts relevant in our context is on Ehrgott (2005); Engau (2007); Ziegler (1995).
A cone in R K is a subset C ⊆ R K such that λ u ∈ C for all u ∈ C and for all λ ∈ R, λ > 0. A cone C ∈ R K is called pointed if u ∈ C implies that (−u) ∈ C for all u = 0. Moreover, a cone C ⊆ R K is called a polyhedral cone if there exists a matrix A ∈ R m×K \ {0} such that C = hcone(A) := {y ∈ R K : A y ≧ 0}. The rows of the matrix A are normal vectors of hyperplanes, and thus a polyhedral cone can be seen as the finite intersection of m (closed and linear) halfspaces. Polyhedral cones can also be described by their extreme rays. This property is an immediate consequence of the well-known Weyl-Minkowski Theorem: Theorem 13 (Weyl-Minkowski-Theorem, cf. Ziegler 1995). A cone C ⊆ R K is finitely generated by n vectors in R K , i.e., C = vcone(B) := Bλ : λ ∈ R n , λ ≧ 0 for some B ∈ R K×n if and only if it is a finite intersection of m halfspaces in R K , i.e., C = hcone(A) = y ∈ R K : A y ≧ 0 for some A ∈ R m×K . Now let C ⊂ R K be a cone. Then the sets are called the dual cone and the strict dual cone of C, respectively. Every cone C ∈ R K induces a binary (ordering) relation R ⊆ R K × R K by defining that (u, v) ∈ R if and only if (v −u) ∈ C. Depending on the context, we also write uRv whenever (u, v) ∈ R (as, for example, in the case of the binary relations defined in (3), (4) and (5) above). Binary relations that are induced by cones are always compatible with scalar multiplication, i.e., (u, v) ∈ R implies (λ u, λ v) ∈ R for all u, v ∈ R K and λ > 0. Moreover, they are compatible with addition, i.e., (u, v) ∈ R implies (u + w, v + w) ∈ R for all u, v, w ∈ R K .
Conversely, binary (ordering) relations that are compatible with scalar multiplication induce cones that represent the respective relation. The following result will be particularly useful in our context.
Lemma 14 (see, e.g., Ehrgott 2005). Let R ⊆ R K × R K be a binary relation on R K which is compatible with scalar multiplication. Then C R := {(v − u) ∈ R K : (u, v) ∈ R} is a cone, and C R induces the binary relation R. If R is additionally compatible with addition, then the following statements hold: It is easy to see that all three orders (3), (4) and (5) are compatible with scalar multiplication and addition. The component-wise order (3) induces the Pareto cone given by P : Similarly, the weak component-wise order (4) induces the closure of the Pareto cone given by R K ≧ = {y ∈ R K : y ≧ 0} = P ∪ {0}, which is a proper, pointed and convex cone (see Lemma 14). Moreover, it is a polyhedral cone that is defined by the identity matrix, i.e., P ∪ {0} = hcone(I) = vcone(I) (where I is the K × K identity matrix). Note also that the Pareto cone is self dual, i.e., P * = P .
Lemma 14 implies that binary relations R that are compatible with scalar multiplications can be equivalently represented by associated cones C R . This interrelation was used, among others, in Engau (2007) to define the concept of cone-efficiency (or C R -efficiency). For a general introduction to ordering cones in the context of vector optimization see, e.g., Tammer & Göpfert (2003); Jahn (2011).
Definition 15 (c.f. Engau 2007). Let Y ⊂ R K be a nonempty set and let C R ⊂ R K be a cone induced by a strict partial order R ⊂ R K × R K (i.e., R is irreflexive and transitive). Then the sets are called the C R -non-dominated and the weakly C R -non-dominated set of Y , respectively. The corresponding pre-images x ∈ X are called C R -efficient and weakly C R -efficient, respectively. Thereby, int(C R ) denotes the interior of C R .
Furthermore, we say that Note that when R is the component-wise order (or Pareto order) defined in (3), then C R = P is the Pareto cone in R K , N (Y, P ) is the non-dominated set (or Pareto set ) of Y , and N w (Y, P ) is the weakly non-dominated set of Y , see, e.g., Ehrgott (2005).

The Ordinal Cone
In the following we show that tail-dominance, represented by the binary relation t , is induced by a polyhedral cone in R K . We will refer to this cone as the ordinal cone. As a first step towards this goal, we prove that Lemma 16. The relation t is compatible with scalar multiplication and with addition.
Proof. To show that t is compatible with scalar multiplication, let λ > 0 and u, v ∈ R K with u t v. It follows that λ Furthermore, it holds that λ u = λ v, and hence λ u t λ v, which implies that t is compatible with scalar multiplication.
It remains to show that t is also compatible with addition. Let u, v, w ∈ R K with u t v, i.e., K i=j u i ≤ K i=j v i for all j = 1, . . . , K and u = v. This implies that K i=j (u i + w i ) ≤ K i=j (v i + w i ) for all j = 1, . . . , K and (u + w) = (v + w), i.e., u + w t v + w. Hence we have proven the compatibility with addition.
It can be proven analogously that ≦ t is also compatible with scalar multiplication and addition.
From Lemma 14 and Lemma 16 we can conclude that the cone C t := {(v −u) ∈ R K : u t v} induced by the strict partial order t is pointed, convex and it does not contain 0. We call this cone the ordinal cone to emphasize that C t equivalently represents ordinal dominance and show that its closure, the cone C t ∪ {0}, is a polyhedral cone that can be described as the intersection of K halfspaces.
Theorem 17. The closure of the ordinal cone is a polyhedral cone. In particular, it holds that C t ∪ {0} = hcone(A t ) with A t ∈ R K×K given by Proof. First note that 0 ∈ (C t ∪ {0}) ∩ hcone(A t ). It thus remains to show that for allũ ∈ R K \ {0}, it holds thatũ ∈ C t if and only if A tũ 0. Now letũ ∈ R K \ {0} with A tũ 0. We define u := 0 ∈ R K and v :=ũ. Hence, it holds v − u =ũ and Thus, we obtain hcone(A t ) \ {0} = C t , which concludes the proof.
Theorem 13 implies that the closure of the ordinal cone C t ∪ {0}, which is a polyhedral cone by Theorem 17, must also have a description based on a finite number of extreme rays. Indeed, the following result provides such a description based on exactly K extreme rays.
Theorem 18. It holds that hcone(A t ) = vcone(B t ) for A t defined according to Theorem 17 and B t ∈ R K×K given by Proof. We first show that hcone(A t ) ⊆ vcone(B t ). Let d ∈ hcone(A t ). Hence, it holds A t d ≧ 0 which is equivalent to K i=j d i ≥ 0 for all j = 1, . . . , K. Set λ j := K i=j d i ≥ 0 and let B j• denote the j-th row of B t , for j = 1, . . . , K.
Note that these descriptions of the closure of the ordinal cone C t ∪ {0} are not unique. Indeed, both the normal vectors in A t as well as the extreme rays in B t could be reordered, and they could be multiplied by arbitrary positive scalars without changing the cone that they define. In the particular description given in Theorems 17 and 18, however, we observe that the matrix B t is the inverse of the matrix A t , i.e., (A t ) −1 = B t .
Remark 19. It holds (hcone(A t )) * = hcone((B t ) ⊤ ) and (vcone(B t )) * = vcone((A t ) ⊤ ) since the normal vectors of the halfspaces given in A t are orthogonal to the extreme rays contained in B t .
It is easy to see that the Pareto cone, P , is a subset of the ordinal cone, C t . Moreover, the dual cone of the ordinal cone is a subset of the Pareto cone, i.e., (C t ) * ⊆ P ⊆ C t . This holds since z ∈ P implies that z 0 and hence A t · z 0, i.e., z ∈ C t . Moreover, z * ∈ (C t ) * is equivalent to (z * ) ⊤ c ≥ 0 for all c ∈ C t which implies z * i ≥ 0, because the i-th unit vector e i ∈ R K is contained in C t for all i = 1, . . . , K. These cones and their duals are visualized in Figure 3.
Remark 20. Note that head-dominance (c.f. equation (1)), represented by the binary relation h , also induces a polyhedral cone. Indeed, we have that Now we can use the ordinal cone C t as described in Definition 15 to reformulate the optimization problem (OCOP) as follows: Here, min C t denotes the minimization in the sense of Definition 15 for the ordinal cone C t = hcone(A t ) \ {0}. In other words, the C t -non-dominated set of problem (OCOP) is given by N (Y, C t ), where Y = c(X). In the following, we use this notation to clearly distinguish between the optimization w.r.t. different ordering cones.

Bijective Linear Transformation Between Ordinal and Pareto Optimization
In the previous subsection we showed that tail-dominance, and hence also ordinal dominance due to Lemma 9, can be equivalently described by the ordinal cone C t . Moreover, the closure C t ∪ {0} of the ordinal cone is the polyhedral cone hcone(A t ) = vcone(B t ) that is spanned by K linearly independent extreme rays in R K , c.f. Theorems 17 and 18. Since the closure of the Pareto cone P ∪ {0} is also a polyhedral cone that is spanned by K linearly independent extreme rays in R K (namely the K unit vectors in R K ), there exists a bijective linear transformation that maps the (closure of the) ordinal cone onto the (closure of the) Pareto cone.
We thus define the following transformed Pareto cone optimization problem where min P denotes the optimization w.r.t. the Pareto cone P according to Definition 15. Note that the objective vector of problem (TOP) corresponds to the incremental tail counting vectorc(x) = A t · c(x) ∈ R K introduced in Section 2.1, that counts in its jth component the number of elements of x that are in category η j or worse. Indeed, for a feasible solution x = {e 1 , . . . , e n } ∈ X we getc(x) = n i=1c (e i ), wherẽ Thus, problem (TOP) is actually a multi-objective optimization problem with K binary objective functionsc 1 , . . . ,c K defined on the ground set S, and with feasible set X ⊆ 2 S . Recall from Section 2.1 thatc 1 (e) = 1 for all e ∈ S and hencec 1 (x) simply counts the number of elements in a solution x ∈ X. Moreover, the vectorc(e) has the consecutive ones property in the sense that whenever a component ofc(e) is zero, then all subsequent components ofc(e) are also zero.
In order to show that the ordinal counting optimization problem (OCOP) (and hence the ordinal optimization problem (OOP)) can be solved by using the above transformation to the "standard" multi-objective optimization problem (TOP), we use a classical non-dominance mapping result for polyehdral cones. This result can be found in Engau (2007) and the references therein among several others. We include a proof, which is similar to the more general proof in Hunt & Wiecek (2003), for the sake of completeness.
Theorem 21 (see, e.g., Engau, 2007). Let Y ⊂ R K be a nonempty set and let hcone(A) be a cone induced by a matrix A ∈ R m×K . Then it holds Proof. This follows immediately from Lemma 9, the relation between orders and cones and Theorem 21.

Solution Strategies
In this section, we discuss a generic algorithmic framework for ordinal optimization problems that takes advantage of the close relationship to multiobjective optimization problems. Since the weighted sum scalarization is a popular approach in multi-objective optimization, the interpretation of weights in the context of ordinal optimization and Pareto optimization is analyzed in more detail, and their relation to numerical representations is discussed.

Ordinal Optimization by Pareto Transformation
From the theory above it follows that we can solve the problems (OCOP) and (OOP) by solving the transformed problem (TOP), which is a standard multi-objective combinatorial problem w.r.t. Pareto optimality. After the computation of the Pareto-efficient set of problem (TOP), or of a minimal complete Pareto-efficient set, respectively, it is necessary to re-compute the corresponding outcome vectors of either problem (OCOP) or (OOP). In this context, a minimal complete Pareto-efficient set of (TOP) is a subset of the Paretoefficient set that contains one Pareto-efficient solution for each Pareto-nondominated outcome vector. We refer to Serafini (1987) for different solution concepts in multi-objective optimization. The efficient sets of the problems (TOP), (OCOP) and (OOP) are equal and will be denoted by X eff in the following. The respective non-dominated sets are denoted by Y TOP nd := N (c(X), P ), Y OCOP nd := N (c(X), C t ) and Y OOP nd , respectively. A procedure for the computation of the efficient set and the respective non-dominated sets based on this Pareto transformation is outlined in Algorithm 1. // ... resp. obj. spaces 5 return efficient set X eff and non-dominated sets Y OCOP nd and Y OOP nd Note that the structural properties of the problems (OCOP) and (OOP) are preserved by the transformation as we do not change the feasible set and as the transformation of the objective function is linear and bijective. In particular, combinatorial solution strategies, like e.g. Bellman's principle of optimality for knapsack problems, can be applied in step 2 of Algorithm 1 to efficiently compute X eff . Ordinal optimization is thus in general no more complex that standard multi-objective optimization.

Weighted Sum Scalarization and Ordinal Weight Space Decomposition
In the following we investigate the interrelation between weighted sum scalarizations for (TOP) and (OCOP) and numerical representations for (OOP). Thereby we rely on the concept of weight space decompositions, which were introduced by Benson & Sun (2000) for multi-objective linear programming and extended to integer linear problems in Przybylski et al. (2010).
The weighted sum scalarization for (TOP) is with λ i > 0 for i = 1, . . . , K and K i=1 λ i = 1. Analogously, the weighted sum scalarization for (OCOP) can be formulated as s is the strict dual cone of the ordinal cone C t , and K i=1 µ i = 1. Recall that the strict dual cone (C t ) * s is the interior of the dual cone (C t ) * , which is visualized in Figures 3(e) and 3(k).
It is a well-known fact that when considering a multi-objective optimization problem, then optimal solutions of weighted sum scalarizations with weighting vectors λ ∈ R K > are always Pareto-efficient (see, e.g., Ehrgott, 2005). Such solutions are called supported efficient solutions. Thus, problem (WSTOP(λ)) always yields Pareto-efficient solutions for problem (TOP). Since (OCOP) can be interpreted as a multi-objective optimization problems w.r.t. the ordering cone C t , every optimal solution of the associated weighted sum problem (WSOCOP(µ)) with weights in the strict dual cone (C t ) * s of C t is ordinally efficient for (OCOP) (see, e.g., Engau, 2007).
The supported efficient solutions of problems (TOP) and (OCOP) are the same, hence there is a one-to-one correspondence between appropriate weighting vectors λ and µ. For a given λ ∈ R K > and x ∈ X we define µ i = i j=1 λ j for all i = 1, . . . , K. Then it holds that which shows that problems (WSOCOP(µ)) and (WSTOP(λ)) have the same objective functions in this case.
Note that µ i = i j=1 λ j and λ i > 0 for all i = 1, . . . , K implies that µ i < µ j for all i < j, as required. Conversely, weighting vectors µ ∈ (C t ) * s satisfy µ i < µ j for all i < j and hence yield associated weighting vectors λ ∈ R K > by setting λ 1 := µ 1 > 0 and λ i := µ i − µ i−1 > 0 for all i = 2, . . . , K. Note also that while the values of µ i = i j=1 λ j , i = 1, . . . , K (for given λ ∈ R K > ) are in general not normalized to satisfy K i=1 µ i = 1, such weighting vectors µ can be easily normalized by setting  The values on dashed lines may not be chosen for λ and µ, because for the weight space decomposition we assume that λ ∈ R K > and 3 i=1 λ i = 1, and for the ordinal weight space decomposition we require 0 < µ 1 < µ 2 < µ 3 and 3 i=1 µ i = 1.
Note that this normalization is applicable since (C t ) * s ⊂ R K > . As a consequence, a weight space decomposition for the multi-objective problem (TOP) can be translated into an associated ordinal weight space decomposition for the ordinal counting optimization problem (OCOP). In this context, a weight space decomposition subdivides the space of relevant weighting vectors λ ∈ R K > with K i=1 λ i = 1 into polyhedral cells such that all weighting vectors from the same cell generate the same efficient solution(s).
Example 24. In the shortest path problem of Example 1 the solutions x 2 , x 3 , x 4 and x 5 are efficient. In Figure 5 the corresponding weight space decomposition is depicted showing the values of λ and µ for which the respective efficient solution is obtained.
From yet another perspective, weighting vectors µ ∈ (C t ) * s , i.e., weighting vectors µ ∈ R K satisfying 0 < µ i < µ j for all i < j, are related to numerical representations as introduced in Section 2.2. Indeed, numerical representations assign a numerical value ν(η i ) to every ordinal category η i , i = 1, . . . , K, such that ν(η i ) < ν(η j ) whenever i < j. Hence we can chose the values µ i , i = 1, . . . , K, equal to the values ν(η i ) of any numerical representation that satisfies ν(η 1 ) > 0. These values can again be normalized without changing the optimal solutions of (WSOCOP(µ)) by setting It is important to note that this does not imply that numerical representations and weighted sum scalarizations are equivalent. Similarly, it is in general not possible to compute all ordinally efficient solutions of problem (OOP) by solving (WSOCOP(µ)) for an appropriate µ. To see this, recall that a solution x ′ ∈ X is called ordinally efficient for problem (OOP) if and only if there is nox ∈ X that ordinally dominates x ′ , i.e., if for everyx ∈ X there exists a numerical representation νx ∈ V such that νx(x ′ ) ≤ νx(x). In contrast, a solution x ′ ∈ X is optimal for problem (WSOCOP(µ)) with appropriate µ if and only if there exists a numerical representation ν * ∈ V such that ν * (x ′ ) ≤ ν * (x) for all x ∈ X.
Remark 25. Ordinal optimization problems may have non-supported efficient solutions. This is illustrated in Example 26. Hence, we can not expect to determine all efficient solutions with the weighted sum method.

Conflicting Real-valued Objectives and Ordinal Objectives
The results of Section 3 can be extended to multi-objective optimization problems that combine a finite number of p "standard" real-valued objective functions w j : X → R with a finite number of r ordinal objective functions o l : X → C l , l = 1, . . . , r, that are in mutual conflict. The number of categories in the l-th ordinal objective function is denoted by K l , i.e., C l = {η l 1 , . . . , η l K l } for l = 1, . . . , r.
For a feasible solution x = {e 1 , . . . , e n } ∈ X, we assume that w j (x) := n i=1 w j (e i ), j = 1, . . . , p. Moreover, o l (x) = sort(o l (e 1 ), . . . , o l (e n )) for l = 1, . . . , r. This leads to the multi-objective ordinal optimization problem with additional cost functions (MOOP): x ∈ X. (MOOP) By replacing the ordered vectors o l (x) for l = 1, . . . , r by the counting vectors c l (x) for l = 1, . . . , r we get a corresponding multi-objective ordinal counting optimization problem with additional cost functions (MCOP): x ∈ X. (MCOP) We denote the concatenated outcome vectors of (MCOP) as v(x) := w 1 (x), . . . , w p (x), (c 1 (x)) ⊤ , . . . , (c r (x)) ⊤ ⊤ ∈ R p+r , wherer := r l=1 K l . Then problem (MCOP) can be transformed into an equivalent standard multi-objective optimization problem w.r.t. Pareto dominance using a linear transformation that is defined by the block diagonal matrix Here, I p×p ∈ R p×p denotes the identity matrix and A l t is the transformation matrix corresponding to the objective c l for l = 1, . . . , r, c.f. Theorem 17. Thus, we get the multi-objective transformed Pareto cone optimization problem (MTOP) Now, problem (MOOP) or, equivalently, problem (MCOP) can be solved by using a simple adaptation of Algorithm 1, c.f. Section 4.

Coherent Real-valued Objectives and Ordinal Objectives
In some practical applications, the elements of S have a real-valued cost (e.g., the length of an edge) and an associated category (e.g., the safety of the corresponding road segment for a cyclist) such that the real-valued cost is, rather than in conflict, coherent with the respective category. This situation is illustrated at the following example: Example 28. Consider the shortest path problem shown in Figure 7. Let w(e) denote the length of an edge e and let o(e) denote its safety: dotted green edges are save and in category η 1 , while solid red edges are insecure and in category η 2 . Then, irrespective of the number of edges contained in the respective paths, the path x 1 = {e 1 , e 2 , e 3 } should be preferred over the path x 2 = {e 4 , e 5 , e 6 } since the total weights are equal w(x 1 ) = w(x 2 ) = 10), and the red sub-path in x 1 has a smaller weight than that of x 2 . In this sense, the weight or length of an edge can be interpreted as an attribute of its respective category. However, x 1 is dominated by x 2 w.r.t. problem (MOOP).
To model the situation where a real-valued objective function w : S → R is in accordance with an ordinal objective function o : S → C with K categories, i.e., the situation where the weight w(e) reflects the multiplicity with which the category o(e) of the element e is to be counted, we introduce a weighted counting The basic idea of this concept is also used in the risk-aware bicycle routing application geovelo, which takes, besides the route length, also the total length of unsafe route segments into account. For example, when K = 4, w(e) = 7 and o(e) = η 2 , then c w (e) = (0, 7, 0, 0) ⊤ . The weighted counting objective of a feasible solution x = {e 1 , . . . , e n } ∈ X equals the sum of the weighted counting vectors of all elements in x, i.e., c w (x) = n i=1 c w (e i ). Thereby, the i-th component of c w (x) corresponds to the total weight of the elements in x that are in category η i , i = 1, . . . , K. Now the weighted counting vector can be handled analogously to the counting vector. Indeed, as in the previous chapter we consider the transformationc w i (x) := K j=i c w j (x) = A t · c w (x) to obtain the weighted transformed Pareto cone optimization problem min Pc w (x) s. t. x ∈ X (WTOP) w.r.t. the concept of Pareto optimality, that can be solved with the methods developed in the preceding sections.

Modelling Aspects
We emphasize that it depends on the context of the respective application whether the multi-objective model (MTOP) or the aggregated model (WTOP) is more suitable. The following example illustrates that the aggregated model (WTOP) is meaningful whenever w and c are interrelated and coherent objectives, while the multi-objective model (MTOP) is particularly useful for unrelated or incompatible objectives.
Example 29. Consider again the shortest path problem depicted in Figure 7. Obviously, the path x 1 = {e 1 , e 2 , e 3 } is the unique efficient solution for problem (WTOP), while the path x 2 = {e 4 , e 5 , e 6 } is the unique efficient solution for problem (MOOP).
Whether x 1 or x 2 are actually preferred thus depends on the interpretation of the weights and of the ordinal categories. Towards this end, suppose that, as in Example 28, the category of an edge corresponds to its security level.
First, consider the case that the real-valued objective w(e) represents the length of the edge e as in Example 28. Then both objectives are interrelated and hence model (WTOP) is appropriate.
If, on the other hand, w(e) represents the toll of the edge or road e, then this real-valued objective is not an attribute of the corresponding category. In other words, both objectives are potentially conflicting and not coherent. Hence, in this case the path x 2 is preferred since the total amount of toll is the same for both paths, but the second path has more green and fewer red edges.

Conclusion
In this paper we investigate ordinal combinatorial optimization problems. We describe different optimality concepts for ordinal objective functions, namely ordinal optimality, which was introduced first by Schäfer et al. (2021), as well as tail-and head-optimality. We prove that all three concepts are equivalent if all feasible solutions have the same length. In general, only ordinal optimality and tail-optimality are equivalent.
We provide alternative descriptions of these three optimality concepts based on associated ordering cones. Using the fact that ordinal optimality and tailoptimality are equivalent, and that tail-optimality can be represented by a polyhedral cone with K extreme rays in R K , we show that ordinal optimization problems can be transformed into equivalent multi-objective optimization problems with binary cost coefficients. The transformation is realized by a bijective linear mapping. The resulting problem can be solved with standard methods from multi-objective optimization, and hence ordinal optimization is as easy or hard as the associated, "standard" multi-objective problems. For example, ordinal knapsack problems and ordinal shortest path problems can be solved by multi-objective dynamic programming, using Bellman's principle of optimality.
The results can be extended to problems with more than one objective function. We suggest two modelling approaches to combine an ordinal objective function with a real-valued objective function. While in the first approach all objectives are considered in a standard multi-objective setting, the second approach allows to model interrelated and coherent objective functions, where the real-valued objective is interpreted as an attribute of the respective category in the ordinal objective.
Future work should focus on the development of tailored optimization algorithms for the associated multi-objective optimization problems that exploit the fact that these problems have binary cost coefficients. Moreover, specific combinatorial problems like, for example, shortest path, knapsack, assignment and general routing and network flow problems should be analyzed both with coherent and with conflicting ordinal and real-valued objective functions.


extfieldsofstudy:
Computer Science
Medicine

text:
Geospatial model of COVID-19 spreading and vaccination with event Gillespie algorithm

We have developed a mathematical model and stochastic numerical simulation for the transmission of COVID-19 and other similar infectious diseases that accounts for the geographic distribution of population density, detailed down to the level of location of individuals, and age-structured contact rates. Our analytical framework includes a surrogate model optimization process to rapidly fit the parameters of the model to the observed epidemic curves for cases, hospitalizations, and deaths. This toolkit (the model, the simulation code, and the optimizer) is a useful tool for policy makers and epidemic response teams, who can use it to forecast epidemic development scenarios in local settings (at the scale of cities to large countries) and design optimal response strategies. The simulation code also enables spatial visualization, where detailed views of epidemic scenarios are displayed directly on maps of population density. The model and simulation also include the vaccination process, which can be tailored to different levels of efficiency and efficacy of different vaccines. We used the developed framework to generate predictions for the spread of COVID-19 in the canton of Geneva, Switzerland, and validated them by comparing the calculated number of cases and recoveries with data from local seroprevalence studies.

Introduction
Numerous epidemic modeling and simulation toolkits have been developed or adapted for COVID-19, ranging from educational models to global-scale comprehensive frameworks [1]. Contact network representations are powerful tools that allow us to understand disease transmission in human populations while taking into account the structure of human interactions, mobility, and contact patterns. Previous studies have paid attention to the influence of various network configurations (scale-free, random, small-world) [2][3][4] on the spreading rate and value of the epidemic threshold. Some models also work with population mobility data, which is especially important in the context of imposing movement restrictions to reduce transmission [1,5].
To improve the accuracy of epidemic models, it is necessary to account for population heterogeneity, e.g., in terms of age, social groups, and mobility patterns, as well as the geographic clustering of infection spread that may result from higher contact rates in locations with higher population density. Examples of models with a detailed representation of these factors include the age-structured household model of Pellis et al. [6] and the two-level (global and local) mixing model proposed by Ball et al. [7] which has further been expanded to account for network structure [8]. However, these models do not include distance metrics or account for differences in population density. In the context of analyzing and predicting the course of the COVID-19 epidemic, the impact of various epidemic control and containment measures needs to be thoroughly examined. [9].

Motivation
The importance of the spatial component in epidemic systems is increasingly recognized. When models are used to evaluate spatially heterogeneous interventions, they must be able to represent the location of hosts and the spatial pattern of transmission in sufficient detail. The inclusion of the spatial component in epidemic models has been considered a complex technical addon that should only be considered when absolutely necessary [10]. However, as software and methods improve and spatial data at the individual level become increasingly available, we can expect the spatial component to play an increasingly important role in infectious disease modeling.

General assumptions
We use a standard representation of the population as a network consisting of N nodes (individuals) {i} and a set of links {e i j }, representing a contact between each node i and j. We assume that the number of individuals in the network remains constant during the simulation [2]. In our case, the network is weighted (all links e i j have a weight depending on the spatial distance between the nodes i and j) and complete (i.e., a link e i j exists for any pair of nodes i and j).
We represent the infection process in the network with a stochastic SEIR model, so that all nodes are at any time point in one of the four states: S (susceptible), E (exposed in the latent period), I (infectious), or R (recovered/removed).
The stochastic infection process S → E occurs as a result of contact between a susceptible (i ∈ S) and an infectious ( j ∈ I ) individual. Two other random processes occur in the network in parallel: the transition of individuals from the exposed to the infectious state E → I and the removal process I → R.
The model assumes two types of contact events, which we call "local" and "global" contacts and distinguish on the basis of the Euclidean distance m i j , with some radius r as the cutoff distance (which is the same for all individuals). Thus, we distinguish between contacts that take place near an individual's residence, and contacts that take place during travel over longer distances.
Local contacts are determined by a cutoff r on the distance m i j . This cutoff is the same for all nodes in the network. Global spread is modeled as a separate process in which each infectious individual can spontaneously transmit the infection to a randomly selected target from the entire population. The process is automatically normalized for population density, as the target is more likely to be located in densely populated areas.
In our model, the following parameters affect the contact probability between two nodes: 1. The assignment of nodes to one of 16 age groups.
The frequency of contacts between age groups is given as an asymmetric 16x16 matrix A = {a i j } [Appendix A]. 2. The distance between nodes. Since only relative distances matter in our model, we use the Euclidean distance between locations on the map. We set a fixed cutoff distance to distinguish between local and global contacts. 3. The population density in a given area. The contact probability for global contacts depends directly on the population density in the area of the contacting nodes-this is achieved by rejection sampling of global contacts over the entire simulated region.

Epidemic model
We represent the stochastic SEIR model as a continuous time three-dimensional Markov chain X = {(S(t), E(t), I (t)) : t ≥ 0} that tracks the number of susceptible, exposed and infectious individuals at each time point. The number of removed individuals can be calculated as N is the total number of individuals, assumed to be constant during the simulation. The epidemic starts with n 0 infectious individuals; the rest of the population is assumed to consist of fully susceptible individuals, i.e., the initial state of X is (S(0), E(0), I (0)) = (N − n 0 , 0, n 0 ). At each time point, the state space of X is described by changing the state of the individuals according to the rules listed in Table 1. All other types of contacts do not change the state of the system.
We use Poisson processes for link activation. When an infection event occurs at a node i ∈ I , we first randomly select a node j within Euclidean distance m i j smaller than r , and activate the link e i j with probability that depends only on the age-contact structure of the population: p i j = a i j /K , where K is the normalization factor: the total contact rate is summed for all age groups. If j is in the susceptible state, the disease is transmitted with probability β L i j = βp i j , where β is the transmission rate and the susceptible node becomes infected (exposed).
For global contacts, we select a random node j over the entire simulated geographic region with rejection sampling, such that the target node is selected according to the population density distribution. This takes into account all non-local contacts (chance encounters, travel, etc.) and allows us to specify their relative proportion with a single parameter: the frequency of global contacts relative to local contacts. In this case, the disease is transmitted with probability β G i j = βp i j d j /D, where d j is the density of the local area of global contact j and D is a normalizing constant.
An infected node remains in the exposed state for a time interval (exponentially distributed) calculated from the mean rate , after which it transitions to an infectious state where it can infect other nodes. The transition to the recovered state follows another Poisson process with exponentially distributed time intervals calculated from the rate γ . A recovered node neither infects other nodes nor is infected by other nodes.

Simulation algorithm
To simulate the infection process, we developed and applied a new variant of the event-modulated Gillespie algorithm. Our implementation supports multiple epidemic control regimes and arbitrary functional forms and distributions of epidemic parameters. It can be used to validate theoretical models and is particularly suitable for simulating epidemic dynamics in large regions such as entire countries.
The unmodified Gillespie algorithm is described as follows. Consider N Poisson processes with rates λ i , i ∈ [1, N ] running in parallel. We denote the density of the event rate for the ith process by ρ i (λ i ). The renewal process is completely characterized by the probability density of inter-event times. For the Poisson process with probability density of event rate ρ i (λ i ), we have In the unmodified Gillespie algorithm applied to epidemic spreading, all possible state transitions for nodes (S → E, E → I , I → R) are Poisson processes with rates i. Multiple simultaneous Poisson processes can be combined into a single global Poisson process with the rate Λ = N i=1 λi, where the time intervals between events are distributed exponentially with the geometric mean 1/Λ, and for each event of this process, we need to calculate what exactly that event is, and where it occurred. (To do this, we need to keep track of all susceptible nodes that can be infected immediately, all exposed nodes that can transition to the infectious state, and all infectious nodes that can be recovered and sample the probabilities of these transitions accordingly.) The original Gillespie algorithm was first developed to track chemical reactions and later adapted to model epidemic processes in well-mixed populations where the only important output is the number of entities in each state (e.g., molecules in chemical models, individuals in epidemic models). However, in network and geospatial models, it is important to know which node exactly makes the transition, leading to computational difficulties.
In our event-modulated Gillespie algorithm, we cut off all unproductive (not leading to a state change) activations, and use a global event queue (a priority queue, i.e., all events are always sorted by time) to automatically place all actual state transitions at the appropriate time points in the future. Therefore, we no longer need to merge all event rates i into the global rate Λ but instead work directly with the SEIR parameters β, , and γ , using them as the rates of the corresponding Poisson processes and generating exponentially distributed time intervals between events.
When a node becomes infectious (either from the exposed state after the incubation period or because it is one of the original infected patients), the infection sequence begins. First, we sample the time interval until removal (recovery or death) as − ln(U )/γ , where U is a random value sampled from the uniform distribution [0.1], and place a removal event in the priority queue. We then generate exposure events and their time offsets, which continue as long as the sum of the time intervals between these events is less than the time between the occurrence of the infection and the removal event. The exposure process is described as follows.
We first determine whether the exposure event will be local or global, the probability of a global transmission being q, and the probability of a local transmission 1 − q. We then sample an age group for the target node from the discrete distribution taken from a row of the contact matrix corresponding to the source node's age group. If the exposed contact is local, we select randomly a node j that is located within the Euclidean distance m ≤ r from the index node i and has an age group equal to the age group we sampled in the previous step. If the exposed contact is global, we select randomly any node j over the entire simulated geographical region using rejection sampling (by generating random uniformly sampled Euclidean coordinates with 10m precision until we have a location that contains a node from the chosen age group), so that the target node is selected according to the population density distribution. Then, we sample the time offset for this exposure event as −ln(U )/β. If the time of the exposure is still less than the time of the recovery, we place the event in the queue and move on to the next exposure event . If the exposure event is already beyond the removal time, we cancel the last exposure event and stop the process.
The continuous time evolution of stochastic processes with known transition rates (including epidemiological applications) can be numerically simulated with the family of Gillespie algorithms, which are statistically accurate. For our needs, we used a particu-lar variation: the event-modulated Gillespie algorithm [11].
The algorithm itself is outlined in (Algorithm 1), adapted from István Z. Kiss et al. [12]. Its core capabilities (infection spreading depending on contact matrices and geographical location of individuals) are abstracted into the find_contact function (line 17). Vaccinations have been modeled as a separate process with piecewise linear intensity, running in parallel with epidemic processes. Vaccination capacities were taken from the vaccination statistics dataset for Geneva provided by the Swiss Federal Office of Public Health FOPH and set at 1000, 5000, and 8000 daily vaccinations (from 31 December 2020, 1 May 2021, and 1 June 2021) as a proxy to represent actual vaccine deliveries and administration. Multiple doses of vaccine are supported (set to 2), with a set interval between doses (set to 21 days, as for Pfizer-BioNTech's COVID-19 vaccine).
Since the simulator works with explicit individuals placed in the 2D environment, we assign a susceptibility value to each individual, initially set to 1.0. Susceptibility functions like a coefficient applied to the probability of infection of the target individual (or equivalently, as a multiplier for the time between infections). Each dose of vaccine administered multiplies the susceptibility coefficient for the target individual by a predetermined number. (Currently set at 0.7-so after a single dose, the individual is 70% less susceptible, and after a second dose, 91% less susceptible, compared to baseline-this represents the efficiency of Pfizer-BioNTech's vaccine. ) We do not explicitly model the difference between vaccine efficiency and vaccine effectiveness-we assume that there is sterilizing immunity and that the spread of virus is stopped if a vaccinated individual has not become infected. (The probability of this is 91% lower after two vaccine doses.) In any case, we consider only the last effect: the reduction in hospitalizations, and fit the simulation accordingly. Currently, the lower efficacy of the vaccines for the new variants is proxied by a correspondingly higher infectiousness.
Simulation supports limiting vaccine availability to selected age groups only. In our simulations, we made vaccines available to people aged 65 and older (starting 31 December 2020), 45 and older (starting 1 May 2021), and 18 and older (starting 1 June 2021).

Simulation running
The main simulation code is written in C++ (C++17) and uses the GDAL [13] library for geospatial transformations and manipulations, and the nanoflann [14] k-d tree implementation for radius queries. The population density map from the Global Human Settlement Layer dataset is loaded as a 2D array, with areas inside and outside the border marked separately. The state of the simulation is exported to an output log on each day of the simulation as current SEIR counts, and as a PNG file representing the current state of the population on the density map, where susceptible individuals are represented by grey pixels, exposed by orange pixels, infected by red pixels, and removed by green pixels. Localized clusters of infection in the simulation as well as the direction of infection spread are thus clearly visible.
A single simulation run for the population of the Canton of Geneva (N=499,480) takes on average 35 minutes. For N=10,000,000, the average run time for a single run was about 120 minutes. The simulation has a high degree of stochasticity, so the time required to complete it varies considerably. A single simulation run consumes about 6 GB RAM for N=10,000,000. All simulation runs are independent and therefore trivially parallelizable. We performed all simulations using the Yggdrasil HPC cluster at the University of Geneva, managed with SLURM.
We run a series of 500 simulations, each simulation starting with the same initial condition where a set of randomly selected nodes is infected and all other N − n 0 nodes are susceptible (in this particular configuration n 0 = 2). We measure the percentage of recovered nodes at the end of the simulation (final size) and use the mean estimator from Seaborn [15] with ±1σ interval to conveniently display the data from multiple simulations. The number of nodes in different states can be aggregated and converted into different indicators to be displayed as stochastic epidemic curves. Figure 1 shows the internal state of the simulation in the population density map, where each pixel corresponds to a single individual, and the epidemic states are shown in different colors.
The distance matrix M is constructed using information about population density for each simulated region, obtained from the ESM2015 dataset [16]. The age-dependent contact matrix A comes from the work of Kiesha Prem et al. [17], which extended the results of the POLYMOD project [18] to 152 countries.

Surrogate model optimization
To provide the initial epidemic parameters for the full simulation, we use surrogate model optimization. We fit a stochastic SEIR model that assumes a homogeneous contact network to match observed data on COVID-19 cases in Switzerland provided by ETH Zurich. [19].
We wrote an adapter function to run the surrogate model and export the simulation results, and compared the results to the actual observed number of hospitalizations for the selected geographic region. Then, we wrote a nonlinear least squares optimization routine curve_fit from the scikit-optimize Python package, which works with arbitrary model functions. We use the trusted region reflected optimization algorithm, which is one of the available algorithms that is general and robust enough to be used in bounded parameters search.
After fitting, the obtained parameter values are used as inputs for the simulation of the main model, leaving only the local infection distance threshold r and the fraction of global contacts in all contacts as free parameters.
In the first stage of the fitting procedure, we fix the values and γ of the SEIR model (inverse incubation time and inverse infectiousness interval, since they do not change) as X and Y , respectively, and assume that β ( contagion rate) changes gradually over time due to the different epidemic control regimes. We assumed that the first infections occurred between January 15 and February 25. The first case was officially registered in the canton of Geneva on 25 February, but the infection had probably already been spreading at that time.
The following parameters were fitted: 1. β k -infection spreading rates for each time interval (β 0 bounded from initial estimations of R 0 ); 2. t 0 -the start date of the epidemics The bounds are shown in Table 2. The fitting routine consists of solving the optimization problem of finding the minimum of the function: ..., θ r ) is a set of parameters we want to estimate, N is the number of available data points, f i (θ ) is the i-th component of the vector of residuals. (In this particular simulation, we are fitting for the number of hospitalizations since this is the most objective and readily available statistic, assuming that the proportion of hospitalized patients always remains the same for each age group, at least within the region under consideration.) Given a model function m(t; θ) and some data points D = {(t I , D I |i = 1, .., N }, we defined the vector of We compared the simulation results with real observations from the clinical dataset COVID-19 provided by the FOPH [19]. The epidemic parameters determined in the fitting step were passed to the main simulation code, which was then run 500 times independently. The initial mean contact radius r for the local transmission mode was set to 500 meters, as adapted from [20], and modified over time between 50 m and 500 m according to changes in mobility. The proportions of global contacts were taken from the same study and ranged from 0.01% to 1.00%, according to the mobility data. For comparison purposes, we used the number of hospitalizations as a reference value. We considered hospitalizations as a subclass of the R (removed) compartment because hospitalized patients do not typically infect other individuals. The percentages of infected persons in each age group assumed to be hospitalized are shown in Fig. 3. We also compared the cumulative percentage of recovered patients in the simulation with figures from seroprevalence studies conducted in the canton of Geneva [21].
The epidemic parameters obtained in the surrogate model are shown in Table 2. The simulation reproduced the first epidemic wave of COVID-19 that occurred in March-April 2020 (Figs. 2 and 3). The simulated number of recovered cases is consistent with seroprevalence measurements conducted in the canton of Geneva in April-May 2020 and November-December 2020 (Fig.  4).

Results
The set of epidemic parameters obtained in the surrogate model is shown in Table 2. The simulation reproduced the initial epidemic wave of COVID-19 that happened in March-April 2020, as well as subsequent periods of rising numbers of cases and hospitalizations (Figs. 1, 2 and 3). Even with the vaccination program in place, the simulation predicts an increase in number of cases in July-August 2021, as new variants of concern are up to 100% more infectious compared to the baseline variant [22]. The effect of vaccinations, however, accumulates steadily leading to the predicted fall in the number of new hospitalizations in mid-September 2021. At the time of this writing, these predictions need to be validated by further observations. Importantly, the simulated numbers of recovered cases is consistent with seroprevalence measurements performed in the canton of Geneva in April-May 2020 and November-December 2020 [21] (Fig. 4).

Results and discussion
We have developed a flexible network simulation model for COVID-19 and successfully tested it with data from the canton of Geneva, Switzerland. Our system allows us not only to track and predict the course of the epidemic, but also to simulate and evaluate different measures put in place to contain the spread of the infection, such as restricting mobility and partial or complete lockdown. The simulation software also makes it possible to model the effectiveness of vaccination and decide which populations (based on age or location) should be prioritized for vaccination.
Because organizations responding to epidemics often have limited resources, optimal resource allocation is critical to develop the most effective interventions. This includes identifying the most vulnerable populations groups, identifying epidemic hotspots, and determining the optimal timing and locations for introducing non-pharmaceutical interventions. Models without a spatial dimension are limited in their ability to answer some questions about intervention effectiveness and resource allocation [8,10]. Our framework Fig. 2 Simulation: daily hospitalizations in the canton of Geneva colored by age group Fig. 3 Daily hospitalizations in the canton of Geneva by age groups (the colored area chart), compared to the simulation results (the black curve) Fig. 4 Percentage of recovered population in the canton of Geneva, aggregated data from 100 simulation experiments with ±1σ error bands, compared to measured seroprevalence estimates. [21] uses age-based contact matrices to account for differences in social networks and contact patterns across age groups. It uses detailed population density maps to predict and visualize geographic locations of epidemic hotspots, and it uses a robust optimization framework to adapt to real-world observations and medical statistics available for the target setting, allowing operators to choose optimal response options.
Our approach also has some limitations. The full representation of the population at the level of individuals is computationally much more expensive than traditional compartmental models. Although the computational resources required are still modest compared to many computational physics problems, this makes continuous experimentation and parameter optimization time-consuming and often impractical. Our approach requires some non-standard parameters that are difficult to calculate from observed data (such as the fraction of "locally" transmitted infections).
The addition of some new features could further improve the simulation and the underlying model. These include accounting for seasonality, dividing the geographic representation into zones (e.g., commercial, residential, business, or recreational areas with corresponding differences in mobility patterns), or accounting for travel across the outer boundaries of the simulated environment (which could be represented as spontaneously generated new infection events). However, the addition of these features risks taking the model far beyond analytical tractability and reducing it to an empirical numerical simulation, limiting our ability to generate new theoretical results. Further research is needed to prioritize the approaches with the greatest robustness and practical utility.
Although this model is more computationally intensive than conventional compartmental models and requires some unconventional parameters, our approach can make the simulation results more reliable. This model provides an efficient way to account for age and geographic structure and mobility without resorting to even more computationally intensive direct simulations of people movement.

Appendix A
See Table 3.


extfieldsofstudy:
Computer Science
Mathematics

text:
FL_PyTorch: optimization research simulator for federated learning

Federated Learning (FL) has emerged as a promising technique for edge devices to collaboratively learn a shared machine learning model while keeping training data locally on the device, thereby removing the need to store and access the full data in the cloud. However, FL is difficult to implement, test and deploy in practice considering heterogeneity in common edge device settings, making it fundamentally hard for researchers to efficiently prototype and test their optimization algorithms. In this work, our aim is to alleviate this problem by introducing FL_PyTorch : a suite of open-source software written in python that builds on top of one the most popular research Deep Learning (DL) framework PyTorch. We built FL_PyTorch as a research simulator for FL to enable fast development, prototyping and experimenting with new and existing FL optimization algorithms. Our system supports abstractions that provide researchers with a sufficient level of flexibility to experiment with existing and novel approaches to advance the state-of-the-art. Furthermore, FL_PyTorch is a simple to use console system, allows to run several clients simultaneously using local CPUs or GPU(s), and even remote compute devices without the need for any distributed implementation provided by the user. FL_PyTorch also offers a Graphical User Interface. For new methods, researchers only provide the centralized implementation of their algorithm. To showcase the possibilities and usefulness of our system, we experiment with several well-known state-of-the-art FL algorithms and a few of the most common FL datasets.

INTRODUCTION
Over the past few years, the continual development of DL capabilities has revolutionised the way we interact with everyday devices. Much of this success depends on the availability of largescale training infrastructures such as large GPU clusters and the ever-increasing demand for vast amounts of training data. In contrary to this, users and providers are becoming increasingly aware of their privacy leakage because of this centralised data collection, leading to the creation of various privacy-preserving initiatives by industry service providers [2] or government regulators [10].
Recently, a viable solution that has the potential to address the aforementioned issue is Federated Learning (FL). FL a term was initially proposed in [26] as an approach to solving learning tasks by a loose federation of mobile devices. However, the underlying concept of training models without data being available in a single location is applicable beyond the originally considered scenario and it turns out to be useful in other practical use-cases. For example, learning from institutional data silos such as hospitals or banks which cannot share data due to confidentiality or legal constraints, or applications in edge devices [18,42].
The main goal of FL is to provide strong privacy protection which is obtained by storing data locally rather than transferring it to the central storage. To solve the underlying machine learning objective each client provides focused updates intended for immediate aggregation. Stronger privacy properties are possible when FL is combined with other technologies such as differential privacy [9] and secure multiparty computation (SMPC) protocols such as secure aggregation [4].
To enable research in federated learning, several frameworks have been proposed including LEAF [7], FedML [15], Flower [5]), (PySyft [37], TensorFlow-Federated (TFF) [19], FATE [42], Clara [29], PaddleFL [25], Open FL [20]. These frameworks are mainly built with a focus to be deployed on real world systems while also providing user with an ability to run experiments with the same code. This desired property often comes with the price that the entry bar for researchers to extend or experiment with these frameworks is limited in a sense that they either need to have extensive experience with distributed systems or require assistance from experts on given framework.
In our work, we decided to take one step back and focus on the construction of a framework that, although not aimed for being deployed to edge devices as primary goal, can provide a useful simulation environment for researchers with the following goals: • Low Entry Bar/ Simplicity: We aim our tool to be as simple as possible for the user while providing all necessary functionalities.
• Extensibility: It is easy to bring your own algorithm or dataset or extend existing ones. We aim to achieve this by providing universal abstractions with a sufficient level of flexibility to experiment with existing and novel approaches to advance the state-of-the-art. • Hardware Utilisation: Cross-device FL experiments are usually of much smaller scale comparing to the centralized setting. This is mainly because of limited device capabilities. Running such experiments on GPU can lead to the under-utilisation of available hardware. We aim to resolve this via the ability to parallelise clients' computation. • Easy Debugging: Debugging multi-process or multi-thread systems is hard. We only require user to provide a single thread implementation which is automatically adjusted to multi-GPU and multi-node setup.
To the best of our knowledge, there is no such tool that could provide a sufficient level of freedom and is simple to use and therefore, we design FL_PyTorch to achieve the aforementioned goal.
FL_PyTorch is an optimization research simulator for FL implemented in Python based on PyTorch [32]. FL_PyTorch is a simple to use tool with its own Graphical User Interface (GUI) implemented in PyQt [33]. During the simulation process, the selected local CPUs/GPUs are accessed in a parallel way. In addition, there is a possibility to use remote CPUs/GPUs. Remote devices are required to have a TCP transport layer for communicating with the master. Regarding supported devices, we target server and desktop stations running on Linux, macOS, or Windows OS for efficient simulations.
For the paper organization, we introduce the general FL minimization objective in Section 2. Subsequently, we provide a deep dive into our FL_PyTorch system in Section 3 and, finally, we demonstrate FL_PyTorch capabilities by concluding several experiments on multiple FL baselines and the most used FL dataset in Section 4.

FL OBJECTIVE AND FEDAVG
In this section, we introduce the FL objective in its general form The global objective is an expectation over local objectives over the randomness inherited from the client distribution P, and the local objectives { } have the form of an expectation over the local datasets {D }. The { } has an additional regularization term that is useful for incorporating prior knowledge of parameters to find .
In the case of the finite number of clients and local data points, both global objectives and local losses { } can be written as simple weighted averages in the empirical risk minimization form (ERM). The most common algorithm to solve (1) is federated averaging FedAVG [26]. This algorithm divides the training process into communication rounds. At the beginning of the -th round ( ≥ 0), the server broadcasts the current global model ( ) to a random subset of clients S ( ) (often uniformly sampled without replacement in simulations). Then, each sampled client performs local SGD updates on its own local dataset and sends its local model update Δ ( ) = ( , ) − ( ) to the server. Finally, the server uses the aggregated model updates to obtain the new global model as follows: where is the relative weight of client . The above procedure will repeat until the algorithm converges. In the cross-silo setting where all clients participate in the training at every round, we have S ( ) = {1,2, . . . , }.

FL OPTIMIZATION SIMULATOR
FL_PyTorch is a system that has been built using Python programming language, and it is based on the DL framework PyTorch. We have made the code repository publicly available 1 . Its backbone consists of a general form of FedAVG displayed in Algorithm 1 partially inspired by Algorithm 1 in [35]. Our proposed general form preserves the standard structure of federated optimization where in each round , subset S ( ) of all available clients is selected. As a next step, the master generates its state which is broadcasted together with the current copy of the global solution ( ) to the selected subset of clients S ( ) . Afterwards, each participating client initializes its local solution ( ,0) to be a copy of the global solution ( ) . Each of these clients then performs steps using the local optimizer ClientOpt with gradient estimated by the LocalGradient function. After this step, each client computes the model and state update, which are then sent back to the master which estimates the global update direction using the ServerGradient function. This estimate is used in ServerOpt that updates the global solution to its new value ( +1) . Lastly, the server global state is updated. In Section 3.2, we show that this general scheme captures all standard algorithms, thus our scheme is sufficient and we believe that it gives researchers a sufficient level of freedom to develop and experiment with each component of our general scheme to push both practical and theoretical FL state-of-the-art.
Our current implementation is fully determined and allows us to configure up to 51 parameters which can be specified either through our easy to use GUI tool or directly via the command line that we discuss in the next subsection. These parameters can be grouped into 4 categories based on their function: • Server Optimizer: the number of communication rounds ; the number of sampled clients per round; global initial local learning rate, global learning rate schedule, and global optimizer with its parameters such as momentum, the name of the algorithm to execute. • Local Optimizer: the number of local epochs or local iterations ; batch size for data loading; local optimizer with its parameters. • Model and Data: model's and dataset's names.
• System Setup: directory to store run metadata, target compute devices, usage of remote compute devices, logging level, number of workers for loading dataset, random seed, thread pool sizes, experiments grouping, user defined comment, enable usage of NVIDIA Ampere GPU Tensor Cores, optional cleanup of PyTorch GPU cache at the end of each round, device for store client state 2 .

Frontend
As we mention previously, FL_PyTorch supports two modes for running federated optimization-Graphical User Interface and Command Line Interface. Below, we discuss both of these in detail.

Graphical User Interface (GUI).
We implemented our Graphical User Interface using the PyQt5 GUI framework [33]. This GUI framework supports all the standard desktop operating systems such as Windows, Linux, or macOS. In addition, the GUI part of FL_PyTorch has built-in VNC server support. If the target machine has not a native windows manager system one can use this mode and connect to GUI part via VNC Viewer software.
To get a better picture of our GUI, we display its main menu provided in Figure 1,  Some extra configurable displayed information from GUI includes a number of currently executing experiments, showing the available size of physical DRAM and swap memory in OS.
Some optimization algorithms require a client to store in its internal state-specific information about algorithm-specific shifts that have a notion about gradient estimator. The analysis is very often irrelevant for the specific strategy of initializing shifts, but in practice, different initialization policies may have consequences in convergence speed. We provide two policies with initializing shift by zero or by full gradient at (0) .
As an example, we show how Algorithm 1 is adjusted to FedAVG and SCAFFOLD. For FedAVG, both global state is an empty dictionary. Locally, we run iteration, which is usually set to be a constant for the theoretical and size of the local dataset for the experiments. LocalGradient returns the unbiased stochastic gradient at ( , ) and ClientOpt is Stochastic Gradient Descent (SGD) with given step-size . Similarly to the global state, local step update is none. ServerGradient is a (weighted) average of local updates and ServerOpt is SGD with fixed step-size 1. For the SCAFFOLD, there are 2 changes comparing to FedAVG. Firstly, the global state ( ) is a non-empty vector of the same dimension as ( ) initialized as zero and each ClientOpt has its own local ( ) , which is set to be a stochastic local gradient at ( ) . LocalGradient then returns the unbiased stochastic gradient at ( , ) shifted by ( ) − ( ) . Secondly, the local model update is sent to master together with the local state update, which is set to the difference between current and previous communication round ( ) = ( ) − ( −1) . The ServerGlobalState is then updated using the average of ( ) 's.

Supported Compressors
In federated learning (especially in the cross-device setting), clients may experience severe network latency and bandwidth limitations. Therefore, practical FL algorithms generally use a communication reduction mechanism to speed up the training procedure. Three common methods to reduce the communication cost are: (1) to reduce the communication frequency by allowing local updates; (2) to reduce communication traffic at the server by limiting the participating clients per round; (3) to reduce communication volume by compressing messages.
We naturally support the first and the second option and we also added support for the compression.
FL_PyTorch allows compressing messages both from the server to local clients and vice-versa. We support several unbiased and biased compressors. The supported unbiased compressor are: Identical compressor (no compression), Lazy or Bernoulli compressor (update is communicated with probability rescaled for the update to be unbiased), Rand-(only random coordinates are preserved uniformly at random rescaled for the update to be unbiased), Natural compressor [17], Standard dithering [1], Natural Dithering [17], TernGrad [41] and QSGD compressor [1]. The supported biased compressor are: Top- [6] and Rank- [38]. For Top-and Rank-compressors, parameter can be specified as a percentage in [0,100] from or as an absolute number greater than 1. Also, we provide means to construct new compressors via function composition and probabilistic switching from existing.   Besides DL models, we provide support for simple quadratic loss function to enable users to explore algorithms in this simplistic regime or debug their implementation. We synthetically generate our objective such that local loss is of the form

Supported models(or pattern) structures
Users are able to specify: dimensionality , strong convexity parameter and smoothness constant for (3) and the number of samples per each client. In addition, we consider two settings wherein the first case all clients optimizes the same objective, i.e., iid setting. In the second scenario, each client has a different objective to optimise, i.e., non-iid setting.

Tracking metrics and experiments analysis
An important feature to evaluate experiments is the ability to track all important quantities during the run of the algorithms. FL_PyTorch allows users to track the number of communication rounds, loss, accuracy, the norm of computed gradient, the norm of the full objective gradient, function values, number of gradient oracle calls, used GPU memory, number of bits which would be sent from workers to master in an actual real-world system, number of clients per round, wall-clock execution time. Furthermore, we provide a dedicated visualization tab in our GUI tool which allows user to load and visualize their experiments in an interactive fashion in several plots. You can observe such quantities during experiment execution via two-dimensional plots. Each experiment has a lot of associated information. We provide means to configure the view of experimental results via selecting needed attributes, sorting experiments by some attribute in the analysis tab, and reordering experiments manually.

The internals of the system
To reach one of our goals of making FL_PyTorch computationally effective, we analyzed various aspects of the PyTorch infrastructure, including the place of PyTorch at NVIDIA computation stack, PyTorch initialization, its forward and backward computation, and speed of typical buses used in a local node setup. This analysis led us to the system design, which asynchronously exploits available hardware.
The Algorithm 1 has several independent parts, and instead of sequential execution of these parts in a simulated environment, independent parts can be partition across independent thread-pools of workers. In our implementation, each worker within a thread pool is a separate CPU thread that lives within a Python interpreter process launched in the operation system (OS) if the system's resources allow it. Since PyTorch and therefore FL_PyTorch allow launching of the computations in a GPUs, for providing separate CPU threads ability to work independently, each thread has its separate GPU CUDA stream to submit its computation work for GPU independently. By design, each worker with the same purpose is assigned to its task-specific thread pool with a thread-safe tasks queue. There are three specific types of tasks -deferred execution, process request for finish execution, and wait on the completion of all submitted work.
Finally, we stress that all this is happening "behind the curtain, " and user interaction, including implementation of the new methods, is agnostic to the aforementioned implementation details.

Bringing custom Algorithm/ Model/ Data
As described in Algorithm 1, each algorithm supported by FL_PyTorch requires implementation of InitializeServerState, ClientState, LocalGradient, ClientOpt, LocalState, ServerOpt, ServerGradient, ServerGlobalState. These are standard centralizedlike PyTorch functions that need to be provided to FL_PyTorch in order to run simulations.
One example of adding new algorithm might be SCAFFOLD assuming that we are given implementation of FedAvg in the required Algorithm 1 format. As described in Section 3.2, one needs to set ClientState to return global shift as proposed in [22] and update LocalGradient to account for both global and local shifts. One further needs to define LocalState to return local shift update, see [22] for the detailed formula. The last modification is in ServerGlobalState, which updates global shift using the local shift updates.
For the datasets, we require users to provide it in a form that is loadable via PyTorch's DataLoaders and for the models, the required format is PyTorch's nn.Module.

EXPERIMENTS
In this section, we provide several experiments to showcase potential use cases of FL_PyTorch . In general, FL_PyTorch can serve as a simulator for researchers to test and compare different algorithms, extend them and analyze each algorithmic component in a plug-and-play fashion as described in Section 3 and Algorithm 1.
For instance, one might be interested in how to reduce the communication burden of a given algorithm. In the literature, three approaches have been proposed. For the first approach, one performs local updates to reduce the frequency of the global model updates, thus also reducing communication frequency. By this strategy, communication cost per client model update can be effectively reduced by a number of local steps. However, while often in practice this leads to improvement, it is not always guaranteed, as it was theoretically shown that there is a tradeoff between communication reduction and convergence [8]. Another possibility is to employ compression techniques that reduce the bits transmitted between the clients and the server. This compression operator should be selected carefully as it comes with an extra variance and it might

SCAFFOLD with Compression
In the first example, we used our framework to analyze whether SCAFFOLD can work with compression, i.e., whether both (Δ ( ) , ( ) ) can be compressed, without a significant decrease in performance. This experiment was carried on a synthetically generated quadratic minimization problem 3. We set the dimensionality of the problem to be 20. Both features 's and responses are generated using uniform distribution on the interval [0, 1]. After this step, we update the data matrix such that the objective is L-smooth with = 2 and strongly convex with = 1 using singular value decomposition (SVD). We generate 10 clients and consider a full participation scenario. The number of communications rounds is chosen to be 100. We set the global learning rate to 0.5 and the local learning rate to 1.0. For the compressor, we choose RAND-, with 3 values of : 100%(no compression), 40% and 20%. We provide results of this experiment in Figure 2. One can observe that dropping 60% coordinates at random has minimal effect on the convergence with respect to iterates. Dropping 80% brings a visible slow-down in per iterate convergence, but it is comparable with respect to the number of communicated bits. This experiment demonstrates that it may be worthwhile to consider an extension of SCAFFOLD by adding compression.

Benefits of Local Updates
For the second experiment, we consider a similar setup as for the previous experiment. We set smoothness constant to 5 and strong convexity parameter = 1. As a compressor we select Lazy/Bernoulli compressor described below4 with = 0.8.
We run MARINA and DIANA algorithms with full gradient estimation for 100 rounds with 10 clients and full participation. These methods were not analyzed in the setting with several local iterations that motivates us to consider these methods for this experiment. Looking into Figure 3, we can observe that both MARINA and DIANA do not benefit from extra local updates. We also experimented with different problem condition numbers ( / ), but it seems that local iterations do not speed up the convergence in the strongly convex regime. This might suggest that a naive combination of these algorithms with local updates will not lead to theoretical benefits, too.

Stochastic Updates and Compression
For the third example, we run MARINA and choose the model to be ResNet18 and for the dataset, federated version of CIFAR10 split uniformly at random (u.a.t.) among 100 clients. We run 200 communications rounds and in each round, we randomly select 25 clients. The local learning rate is 1.0 and the global learning rate is 0.001. We perform one local iteration, i.e., = 1.
We investigate the effect of the combined effect of compression and stochastic gradient estimation since both of these approaches introduce extra variance that might negatively impact the quality of the obtained solution.
For gradient estimation, we employ SGD-NICE sampling strategy that estimates gradient via selecting uniformly at random a subset of samples of fixed size.
in the experiments reflects the relative subset size with respect to the total size of the local dataset. The compressor operator is Bernoulli compression. Looking into Figure 4, one can see that having Lazy compression with = 0.1 combined with SGD-NICE with = 10% hurts the convergence. On the other hand, having a less aggressive compressor or using a bigger sample size for SGD-NICE does not deteriorate the performance.

CONCLUSIONS
In this work, we have introduced FL_PyTorch , an efficient FL simulator based on PyTorch that enables FL researchers to experiment with optimization algorithms to advance current state-ofthe-art (SOTA). FL_PyTorch is an easy to use tool that supports SOTA FL algorithms and the most used image classification FL datasets. In addition, FL_PyTorch uses several levels of parallelism for efficient execution while being simple to extend using only standard PyTorch abstractions. In its basic version, it does not require any programming and all its pre-implemented components are accessible through its GUI which allows to set up, run, monitor and evaluate different FL methods.
For the future, we plan to provide users with more pre-implemented algorithms, datasets, and models. We also aim to extend our visualization tool. After open-sourcing FL_PyTorch , we expect a significant increase in the number of supported algorithms, datasets, and models via inputs from the FL research community.  Step-1: Master -Initialize Model -Prepare shared Datasets for each thread in a training Thread Pool -Send request for local optimization handled with a local thread pool Step-3: Master Does not perform waiting of workers completion Step-4: Master -Compute weight for aggregation based on data statistics -Get the next ready model for aggregation -Aggregate Step-4: Master Step-5: Completion