{'fraction_non_alphanumeric': 0.044934377319905315, 'fraction_numerical': 0.011485719888504475, 'mean_word_length': 4.046939064046374, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 37, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Novel Machine Learning Approach to Data Inconsistency with respect to a Fuzzy Relation

Inconsistency in prediction problems occurs when instances that relate in a certain way on condition attributes, do not follow the same relation on the decision attribute. For example, in ordinal classification with monotonicity constraints, it occurs when an instance dominating another instance on condition attributes has been assigned to a worse decision class. It typically appears as a result of perturbation in data caused by incomplete knowledge (missing attributes) or by random effects that occur during data generation (instability in the assessment of decision attribute values). Inconsistencies with respect to a crisp preorder relation (expressing either dominance or indiscernibility between instances) can be handled using symbolic approaches like rough set theory and by using statistical/machine learning approaches that involve optimization methods. Fuzzy rough sets can also be seen as a symbolic approach to inconsistency handling with respect to a fuzzy relation. In this article, we introduce a new machine learning method for inconsistency handling with respect to a fuzzy preorder relation. The novel approach is motivated by the existing machine learning approach used for crisp relations. We provide statistical foundations for it and develop optimization procedures that can be used to eliminate inconsistencies. The article also proves important properties and contains didactic examples of those procedures.


Introduction
Ordinal classification (also called ordinal regression) problems constitute a very important part of machine learning and statistical analysis [1]. In ordinal classification, the goal is to predict for a certain instance u from set U, one of K different ordinal class labels y ∈ {1, . . . , K}. Usually, u is characterised by its values for a given set of condition attributes, while y is called a decision attribute. Ordinal classification problems exploit the existing ordering on the decision attribute. In some cases, an ordering also exists on the condition attributes. One way to incorporate that knowledge is through socalled monotonicity constraints. For a given preorder (dominance) relation on the set of instances U based on the condition attributes, the monotonicity constraints can be formulated as follows: if instance u 1 dominates u 2 w.r.t. the given dominance relation on the condition attributes, then u 1 should be assigned to the same or to a better decision class than u 2 . In such a case, we say that u 1 is consistent with u 2 . Obviously, consistency is a symmetric relation and instances that are incomparable w.r.t. condition attributes are consistent by default.
Ordinal classification problems that include monotonicity constraints are called monotone classification problems. They arise in many areas, such as medical diagnosis [2], survey data [3], bankruptcy risk estimation [4], house pricing [5] and others. A comprehensive survey of monotone classification methods is given in [6].
In practice, not all pairs of instances satisfy the monotonicity constraints (are not consistent) due to some imperfectness of ordinal classification data, like missing attributes or instability of the assessment of decision attribute values at the stage of data generation. Hence, some initial preprocessing of data is often performed in order to enforce the monotonicity constraints. One of the main methods to handle inconsistencies in data is the rough set approach [7]. For a given decision class in a classification problem, the approach outputs lower and upper approximations of that class. The lower approximation contains instances from the decision class that are consistent with all other instances, while the upper approximation contains instances that relate with instances from the decision class. The original rough set approach handles inconsistencies w.r.t. an equivalence (indiscernibility) relation and not w.r.t. a dominance relation (as is defined in the monotonicity constaints). To make it applicable for monotone classification problems, Greco et al. [8] extended the original rough set theory with their Dominance-based Rough Set Approach (DRSA) where the new approach replaces the indiscernibility relation by a dominance relation. After the introduction of DRSA, the original rough set theory is usually referred to as Indiscernibility-based Rough Set Approach (IRSA). Recently, the two were integrated into the Preorder-based Rough Set Approach (PRSA) [9].
A more comprehensive analysis of monotone classification from the statistical learning point of view was given by Kot lowski and S lowiński [10]. They provided statistical foundations of the monotonicity constraints and developed a machine learning method to incorporate them into data analysis. This method removes inconsistencies in data ("monotonizes" them) in result of an optimization procedure that minimizes the cost of label changes in the decision attribute. It produces a new set of labels called a monotone approximation. This approach also generalizes standard rough sets and provides another probabilistic view of them. The approach found its application in the same areas as DRSA [11], as well as in the development of rule induction and ensemble rules methods [12]. In the remainder of the paper, we refer to this method as KS approach.
On the other hand, fuzzy logic and fuzzy set theory [13] study the gradual truth of logical statements, and are used extensively in modeling imprecise and vague information. One of the ways to utilize fuzzy logic in data analysis is through fuzzy relations that are able to model relationships between numerical vectors. Namely, the usual crisp relations may distinguish only between two extreme cases: either instances relate or do not. Fuzzy relations, on the other hand, can express a degree in which two instances relate on a scale between 0 and 1. This is suitable to model similarity between numerical vectors or other structures (graphs, strings, DNA chains ...). The integration of fuzzy logic and IRSA was initially proposed by Dubois and Prade [14], allowing to approximate fuzzy sets using a fuzzy indiscernibility relation. A similar extension of DRSA to fuzzy set theory was proposed by Greco et al. [15]. This article is motivated by the KS approach introduced above in the sense that we generalize the monotonicity constraints using fuzzy relations while the ordinal classes are replaced with fuzzy membership degrees. Instead of a crisp preorder relation (or an equivalence relation if it is symmetric), we will now consider a fuzzy T -preorder relation to model the relationship between different instances on the condition attributes, where T refers to a given t-norm that models conjunction in fuzzy logic. The T -preorder relations also include T -equivalence relations that can measure (symmetric) similarity between numerical vectors. Moreover, the new approach requires that the decision attribute is a fuzzy set, i.e., it has to take values from interval [0, 1]. Hence, it is appropriate for problems where the decision attribute can be modeled using values from this interval; concretely, for binary classification and regression problems.
Just like the KS approach [10], our proposal is also interesting from the granular computing point of view. Granular computing is a paradigm which involves a partition of information into meaningful groups, classes or clusters called granules [16,17,18], and which has been applied to diverse models in data analysis. For example, in [19] and [20], granular computing using neighborhood systems for the interpretation of granules was studied, while in [21,22,23,24], granular aspects of rough set theory were examined.
In particular, the sets obtained with the Kot lowski-S lowiński approach, as well as with the novel approach, possess the property of granular representation: they can be represented as unions of meaningful granules [9,25,26]. Such sets are called granularly representable sets [9]. Due to the granular properties of our new approach, we call its result a granular approximation.
Granularity of fuzzy rough sets has been already applied in rule induction [27], and we expect that the granular approximations we propose can serve the same purpose. Since granular approximations are natural generalizations of fuzzy rough sets, they can be applied in methods like Fuzzy-Rough Nearest Neighbours (FRNN) for classification [28,29], Fuzzy Rough Prototype Selection (FRPS) [30], Fuzzy Rough Feature Selection (FRFS) [31,32], and so on.
The article is also an example of a successful integration of ideas and contributions of rough sets, fuzzy sets and machine learning. Handling inconsistency and granulation are main contributions of rough sets. The theory of fuzzy sets allows us to use fuzzy relations to model a non-binary interaction among instances. For example, fuzzy relations can be used in modeling similarities between instances represented by numerical vectors. At the end, including statistical/machine learning allows us to make data consistent, incurring the least possible cost (w.r.t. some loss function) using optimization methods.
The remainder of the paper is organized as follows. In Section 2, we recall the required preliminaries about statistical learning theory, monotone approximations, fuzzy logic, and fuzzy rough sets. In Section 3, we develop the statistical foundations of granular approximations. Section 4 deals with optimization problems that output granular approximations, while their important properties with proofs are given in Section 5. Section 6 provides didactic examples for the approaches from Section 4. Section 7 contains our conclusion and outlines future work.
In Appendix A-D, we deal with the dual formulations of the optimization problems introduced in Section 4. Using the duality theory, we obtain greedy algorithms for the optimization problems from Section 4 that allow us to prove Proposition 5.2.

Statistical learning for monotone classification
A random variable X is a mapping from a probability space to a certain codomain X. If the codomain is a subset of the real numbers, X is usually characterized with a cumulative distribution function (CDF) defined as F X = P (X ≤ x) for x ∈ X. A CDF is a non-decreasing and right-continuous function with codomain [0, 1]. If the CDF is continuous then we say that X is continuous, while if the image of the CDF is a finite set, we say that X is discrete. Based on the CDF, a quantile function may be defined as follows: Q X (p) = inf{y; F X (y) ≥ p} for 0 < p < 1. In other words, if p is in the image of F X , then Q X (p) is the smallest value for which P (X ≤ Q X (p)) = p. The value Q X ( 1 2 ) is called the median of X . The expected value of X can be expressed using the quantile function [33]: We say that X 1 stochastically dominates Proposition 2.1. [34] For two random variables X 1 and X 2 , it holds that The above proposition states that the stochastic dominance can be characterized using quantile functions instead of CDFs.
We now examine the prediction problem. Let X and Y be two random variables with codomains X and Y respectively. When making predictions, we examine a causal relationship between X and Y, i.e., how does X influence Y. Concretely, we are interested to find a function h such that h(X ) is close to Y, i.e., it predicts values of Y for given values of X . Formally, let L : Y × Y → R + be a loss function. A prediction problem consists in finding a function h : X → Y such that the risk is minimized. The optimal h, denoted as h * , is called the Bayes predictor. In practice, the random variables X and Y are unknown and we only have their realizations x 1 , . . . , x n and y 1 , . . . , y n . Our goal is then to minimize the empirical risk:R Minimization of the empirical risk is called learning and it basically amounts to an estimation of the Bayes predictor. The causal relationship between X and Y may be represented by a family of random variables Y X =x , which stands for variable Y conditioned on X = x. Such a random variable, for a fixed x, may be described by its CDF: Searching for an optimal prediction function h in the learning process may be seen as an estimation of certain characteristics of the family of random variables Y X =x . For example, when the loss function is mean squared error for y,ŷ ∈ Y and Y = R, then the Bayes predictor is h * (x) = E(Y|X = x), i.e., the conditional mean, while if the loss function is mean absolute error L(y,ŷ) = |y −ŷ|, then the Bayes predictor is h * (x) = Q Y|X =x ( 1 2 ), i.e., the conditional median [35]. The previous examples show that a Bayes predictor is a characteristic of family Y X =x (conditional mean and median in the examples), which means that the learning process leads to an estimation of those characteristics.
Kot lowski and S lowiński [10] introduced a statistical framework for monotone classification. In this case, it is assumed that there is a preorder relation X on codomain X of X while Y consists of a finite number of totally ordered values that distinguish different ordinal classes. Denote these classes by 1, . . . , K. The monotonicity constraint states that if x x ′ then x has to belong to at least the same class as x ′ . This is also called the Pareto principle in decision theory. Let K −1 = {1, . . . , K − 1}. In probabilistic terms, the monotonicity constraint says that x x ′ implies The previous expression means that the probability that x will be assigned to a class at most k is smaller or equal than that x ′ will be assigned to the . The goal of monotone classification is to find a proper monotone h under the assumption that the family Y X =x is monotonically constrained. Since h, as the output of the learning process, should be as close as possible to the Bayes predictor h * , we require that h * is also monotone. Given that the form of h * depends on the loss function, choosing a proper loss function is crucial for the learning process. A loss function for which the Bayes predictor is monotone is called a monotone loss function. Kot lowski and S lowiński [10] showed that both mean squared error and mean absolute error are monotone loss functions. They also examined a family of monotone loss functions called p-quantile loss defined as: for p ∈ [0, 1], where 1 stands for the indicator function. The name p-quantile loss is used since the Bayes predictor for such loss function is the conditional p-quantile h * p (x) = Q Y|X =x (p). For p = 1 2 we have that L 1/2 is equivalent to the mean absolute error. For the p-quantile loss, we have the following important result proved in [10].
Proposition 2.2. Let s : Y → R be an increasing function. Then the loss functions L p (y,ŷ) and L p (s(y), s(ŷ)) have the same Bayes predictor. Proposition 2.2 states that a different scaling of ordinal classes does not change the Bayes predictor, only the order matters.

Monotone approximation
In order to incorporate monotonicity constraints into the learning process, KS approach use an optimization procedure to "monotonize" data by eliminating inconsistencies. Letȳ i , i = 1, . . . , n, be the observed ordinal labels which do not necessarily satisfy monotonicity constraints due to possible inconsistency, and letŷ i , i = 1, . . . , n, be the values that we want to learn and which satisfy the constraints. Then, for a given monotone loss function L, the optimization problem can be formulated as In other words, one wants to calculate new labels that are as close as possible to the original ones w.r.t. loss function L and which satisfy the monotonicity constraints. The obtained labels are called a monotone approximation of the original ones. The same authors showed that when L is monotone, then problem (5) can be solved using linear programming. Moreover, the solutions of the linear optimization problem will always be integers due to the unimodularity of the constraint matrix [36].

Fuzzy logic connectives
The definitions and terminology in this subsection are based on [37]. Recall that a t-norm T : [0, 1] 2 → [0, 1] is a binary operator which is commutative, associative, non-decreasing in both arguments, and for which it holds that ∀x ∈ [0, 1], T (x, 1) = x. Since a t-norm is associative, we may extend it unambiguously to a [0, 1] n → [0, 1] mapping for any n > 2. Some commonly used t-norms are listed in Table 1.
We say that x ∈ [0, 1] is a nilpotent element of a t-norm T if there exists a natural number n such that A t-norm is called nilpotent if it is continuous and every x ∈ (0, 1) is a nilpotent element. For example, T L from Table 1 is nilpotent while the others are not. A t-norm is strict if it is continuous and strictly increasing in both arguments. T P from Table 1 is strict while the others are not.
T P , T L and T D from Table 1 are Archimedean, while T M and T nM are not. Table 1: Some common t-norms and their R-implicators A t-norm is a continuous Archimedean t-norm if and only if it is either strict or nilpotent.
We denote with a family of nilpotent t-norms, i.e., t-norms that are isomorphic to T L with bijection ϕ and we denote with a family of strict t-norms, i.e., t-norms that are isomorphic to T L with bijection ϕ.
We say that a t-norm is D-convex if its partial mappings are convex functions. T P and T L from Table 1 are D-convex, while the others are not.
More details on the characterization of D-convex t-norms can be found in [9].
The residuation property holds for a t-norm T and an implicator I if It is satisfied if and only if T is left-continuous and I is defined as the residual implicator (R-implicator) of T , that is The very right column of Table 1 shows the residual implicators of the corresponding t-norms. Note that all of them, except I T D , satisfy the residuation property. Implicators that satisfy the residuation principle have the ordering property x ≤ y ⇔ I(x, y) = 1.
Given a [0, 1] → [0, 1] bijection ϕ, the residual implicators of nilpotent and strict t-norms T L,ϕ and T P,ϕ will be denoted by I L,ϕ and I P,ϕ . For a left continuous t-norm T and its R-implicator I, we define the negator induced by I as N(x) = I(x, 0). We will call a triplet (T, I, N) obtained as previously explained a residual triplet. If a t-norm from a residual triplet is continuous and Archimedean, then the negator of the triplet is involutive if and only if the t-norm is nilpotent. In such case, the negator has the form Proposition 2.4. For a residual triplet, the following holds: I(T (x, y), z) = I(x, I(y, z)).
The standard negator is obtained when, for example, the t-norm is the Lukasiewicz one. In general, a t-norm for which the induced negator of its R-implicator is involutive is called an IMTL t-norm. We will call a residual triplet (T, I, N), that is generated with an IMTL t-norm, an IMTL triplet.

Fuzzy sets and fuzzy relations
Given a non-empty set U, a fuzzy set A in U is an ordered pair (U, m A ), where m A : U → [0, 1] is a membership function that indicates how much an element from U is contained in A. Instead of m A (u), the membership degree is often written as A(u). If the image of m A is {0, 1} then we obtain a crisp or classical set. For a negator N, the fuzzy complement coA is defined as coA(u) = N(A(u)) for u ∈ U. If A is crisp then coA reduces to the standard complement. For α ∈ (0, 1], the α-level set of fuzzy set A is a crisp set defined as A α = {u ∈ U; A(u) ≥ α}.
A fuzzy relation R on U is a fuzzy set on U × U, i.e., a mapping R : U × U → [0, 1] which indicates how much two elements from U are related. Some relevant properties of fuzzy relations include: A reflexive and T -transitive fuzzy relation is called a T -preorder relation while a symmetric T -preorder relation is called a T -equivalence.

Fuzzy rough and granular approximations
Let U be a set of instances, A a fuzzy set on U and let R be a T -preorder relation on U. The fuzzy PRSA lower and upper approximations of A are fuzzy sets for which the membership function is defined as: for u ∈ U. The approximations have some important properties [38]: Given λ ∈ [0, 1], T -preorder R, t-norm T and u ∈ U, a fuzzy granule is defined as a parametric fuzzy set where the union of fuzzy sets is defined with the max operator.

Ontic fuzzy sets and probabilistic uncertainty
Fuzzy sets are often related to uncertainty modeling [39,40,41]. However, we should be very careful when mentioning that fuzzy sets are used to model uncertainty. First, two types of classical (crisp) sets have to be distinguished: conjunctive and disjunctive sets [42]. A conjunctive set is a collection of items that represents a well known complex entity, i.e., it is a conjunction of its elements. For example, a time interval that describes a span of some activity. On the other hand, a disjunctive set describes incomplete information about an ill-known object. The object of interest is contained in the disjunctive set but we do not know which element it is, i.e., the set is a disjunction of its elements. For example, an event that occurred at an unknown moment in time is described with a time interval that represents our knowledge about the unknown event. Conjunctive sets are also known as ontic sets while disjunctive sets are called epistemic sets. Fuzzy sets are used to model gradual information which is not uncertain by itself. Fuzzy sets may be related to uncertainty only if the underlying universe, on which a fuzzy set is defined, is a disjunctive set. In that case, fuzzy sets make incomplete knowledge more expressive by allowing gradual information. Such fuzzy sets are usually known as epistemic fuzzy sets and form the basis of possibility theory [43]. In this article, we always use fuzzy sets defined over a conjunctive universe, i.e., ontic fuzzy sets, while we assume that the uncertainty in data is of probabilistic nature. An example of an ontic fuzzy set is a set of people that are "tall", i.e., a fuzzy set whose universe is some set of humans and its membership function is a height measure of those humans. The height is an actual physical characteristic. In this case, no uncertainty or lack of knowledge exists.

Granularly representable random fuzzy sets
We assume that we observed a finite set of instances U from the underlying universe, i.e., U is a random sample. U is described with condition and decision attributes where the decision attribute takes values in [0,1], which are interpreted as membership degrees to an unknown fuzzy set that we want to reconstruct using the observed values. From the perspective of statistical learning theory introduced in the previous section, condition attributes correspond to random variable X while the decision attribute corresponds to random variable Y, which now takes values from interval [0, 1]. The fuzzy set that we want to reconstruct contains uncertainties that are represented in a probabilistic way, i.e., we assume that the actual values are altered due to perturbation. Perturbation may be caused by the incomplete knowledge about data (missing attributes) or by random effects that occur during data generation. Such altered values are represented by a family of random variables {A(u), u ∈ U} which model our uncertainty about the ill-known membership degrees {A(u), u ∈ U}. In other words, for each instance u, the ill-known membership degree A(u) is represented with the random variable A(u) having codomain [0, 1]. Family {A(u), u ∈ U} is a special case of a random fuzzy set defined in [44] (the other name is fuzzy random variable). Hence, we may refer to the family as random fuzzy set A.
The family {A(u), u ∈ U} corresponds to family Y X =x from the previous section. Therefore, we formulate the reconstruction of fuzzy set A as the problem where for a given set of instances U and its condition and decision attributes, we want to estimate characteristics of A(u) (like conditional mean, median and quantiles mentioned above) in order to describe the illknown A(u). Knowledge about condition attributes is represented using a T -preorder relation R, i.e., for each pair u, v ∈ U we are given the value R(u, v). We denote the observed decision values asĀ(u) for u ∈ U.
In the first step, we will extend the monotonicity constraint (3) for a T -preorder relation. We have the following simple result.
Proposition 3.1. A fuzzy set A in U is granularly representable if and only if it satisfies the object monotonicity property, i.e., The second equivalence holds from the observation that the maximum is reached for u = v due to reflexivity of R. The third equivalence holds from the residuation property. Proposition 3.1 reveals that the granular representability is equivalent to the object monotonicity property. On the other hand, the object monotonicity property is the fuzzy extension of the crisp Pareto principle whose probabilistic form is given in (3). That can be seen from the fact that inequality T ( R(v, u), A(u)) ≤ A(v) is the fuzzy expression of the sentence: "If v dominates u and u ∈ A then v ∈ A". We extend the principle provided in (3) for the case when the preorder (dominance) relation " " is replaced with a fuzzy T -preorder relation R and when the decision attribute is expressed in the form of membership degrees rather than binary class membership. For the terminology, instead of "monotonically constrained" we use the term "granularly representable" introduced in [9]. The new term is more appropriate because a T -preorder relation is more general and may express different things, not only monotonicity (e.g., if the relation is a T -equivalence, then we are talking about similarity). In order to relate the granular representability and the family of random variables {A(u), u ∈ U}, we introduce the following definition.
Definition 3.1 is an extension of the third equivalence in (3). It states that A is granularly representable if all its p-quantiles A p (p ∈ [0, 1]) are granularly representable as ordinal fuzzy sets.
The next question is, if the random fuzzy set A is granularly representable, is its expected value EA, defined as EA = {E(A(u)), u ∈ U}, also granularly representable? Before answering this question, we recall the well-known Jensen inequality [45].
Proposition 3.2. Let µ be a probability measure on the set of reals, g a µ-measurable real function and φ a real convex function. It holds that Since the standard (Lebesgue) measure is equivalent to the probability measure on [0, 1] (measure value of interval [0, 1] is 1), the above inequality translates to Using Jensen's inequality, we obtain the following result.
Proposition 3.3. Let T be a D-convex t-norm and I its R-implicator. Then EA is granularly representable as soon as A is.
Proof. For every u, v ∈ U, we need to prove that Using (1), we have that ∀u ∈ U, EA(u) = 1 0 A p (u)dp. It follows that The first inequality follows from the fact that T (c, ·) is a convex function for a constant c and Jensen's inequality. The second inequality follows from the granularity of A p .

Calculation of granular approximations
In this section, we discuss which properties of A can be estimated and how to do this in practice. In general, the observed fuzzy setĀ is not granularly representable due to the presence of inconsistency, so our goal is to find a granularly representable set that is close to it by minimizing a certain loss function. For given loss function L, the general form of our optimization problem is where {Â(u), u ∈ U} is the unknown granularly representable set. We will call the result of optimization problem (11) the granular approximation of fuzzy set {Ā(u), u ∈ U}.
Optimization problem (11) is the main contribution of the article. It allows us to remove inconsistencies (obtain a granularly representable set) with the least cost of alteration of values (w.r.t. loss function L). The remainder of the section investigates specific cases for which problem (11) can be efficiently solved.
Under the assumption that A is granularly representable, it is desirable to use loss functions for which the Bayes predictor is granularly representable as well.
Definition 4.1. We say that a loss function L is granular with respect to a left-continuous t-norm T and fuzzy relation R if its Bayes predictor is granularly representable under the assumption that the underlying family of random variables {A(u), u ∈ U} is granularly representable w.r.t. T and R.
Note that with this definition, the p-quantile loss function (4) is granular, since its Bayes predictor is the quantile fuzzy set A p , which is granularly representable by definition. The mean squared error loss function (2) is granular for D-convex t-norms since the Bayes predictor EA is granularly representable in this case by Proposition 3.3. Hence, both loss functions introduced in Subsection 2.1 are suitable for the calculation of granular approximations.
In problem (11), both objective function and constraints are not necessarily linear and may take different forms that depend on loss function L and on the type of fuzzy logic connectives used. However, in case of the loss functions (4) and (2), and continuous Archimedean t-norms, the optimization problem can be efficiently solved.
Indeed, consider t-norms T L,ϕ and T P,ϕ introduced in Eq. (6) and (7). If T L,ϕ is used in (11), then the set of constraints that express granular representability can be simplified in the following way: for all u, v ∈ U, where we introduced the shorthands R ϕ (u, v) = ϕ( R(u, v)), α u = ϕ(Â(u)) and α v = ϕ(Â(v)). The last equivalence holds because 0 is always smaller than α u , hence the max operator can be lifted.
If T P,ϕ is used then in an analogous way we find The border constraints now become 0 ≤ α u ≤ 1 for all u ∈ U. We can conclude that using continuous Archimedean t-norms leads to linear optimization constraints. This is a promising result since many optimization solvers are very efficient with linear constraints.
In both cases, the empirical risk can be expressed as In the empirical risk above, the non-linear term ϕ −1 (α u ) appears. Function ϕ −1 is an arbitrary bijection which can lead to a non-convex optimization problem. However, Proposition 2.2 states that a different scaling of values does not change the Bayes predictor delivered by the p-quantile loss function.
To eliminate the non-linearity, we can apply ϕ to both parameters of the loss function and replace L p (Ā(u), ϕ −1 (α u )) by L p (ϕ(Ā(u)), α u ). Although the value of the estimand (the quantity that is estimated, i.e., the Bayes predictor A p ) remains unchanged with the new loss function, the estimator (the result of the optimizationÂ p ) can be different. From the theory of quantile regression, we can express the optimization of the quantile risk as a linear program [46]. We introduce new variables x u , u ∈ U and y u , u ∈ U such that , 0), as well as the shorthandĀ ϕ (u) = ϕ(Ā(u)). In case T L,ϕ is used, we can reformulate optimization problem (11) as In case of T P,ϕ , optimization problem (11) obtains the form Summarizing, for quantile risk and a continuous Archimedean t-norm, the optimization problem (11) can be expressed as a linear program and, therefore, efficiently solved using one of many existing efficient linear programming solvers. We have the following technical result. (12) and (13), are redundant.
Proof. Assume that the constraints are removed and that an optimal solution α * u , u ∈ U, has values smaller than 0 or larger than 1. We construct another solution from α * u , u ∈ U, by replacing values larger than 1 by 1, and values smaller than 0 by 0. It is easy to check that the new solution satisfies the granular representability constraints. From the constraints x u −y u =Ā ϕ (u)− α u , u ∈ U, it is easy to see that when α u ≥ 1 thenĀ ϕ (u) − α u ≤ 0, which leads to x u = 0 and y u = α u −Ā ϕ (u), and when α u ≤ 0 thenĀ ϕ (u) − α u ≥ 0, which leads to x u =Ā ϕ (u) − α u and y u = 0. Hence, after replacing values larger than 1 by 1, the values of y u will be reduced and after replacing values smaller than 0 by 0, the values of x u will also be reduced. In both cases, the value of the objective function will be reduced. Therefore, we constructed a feasible solution with a smaller cost which contradicts the optimality of α * u , u ∈ U.
A solution of optimization problems (12) and (13) is not necessarily unique. The reason is that the estimated CDF of A(u) for some u ∈ U is not necessarily increasing and may have a region where it is constant. In such case, we obtain an uncountable family of optimal solutions. However, if for some probability parameter p we have infinitely many solutions, the lower and upper bounds of such family of solutions can be calculated by running the linear programs with parameters p − ǫ and p + ǫ, respectively, for sufficiently small ǫ.
If the mean squared error is used as a loss function, it is obvious that the objective function will become non-linear. Also, Proposition 2.2 does not hold anymore and using L p (Ā ϕ (u), α u ) instead of L p (Ā(u), ϕ −1 (α u )) will lead to the estimation of a different Bayes predictor. However, we will include this approach in our analysis since it may give good results in practical applications. In this case, the optimization problem for the t-norm T L,ϕ is while for T P,ϕ the corresponding problem is Using a similar argument as in Proposition 4.1, we may drop the constraints 0 ≤ α u ≤ 1, u ∈ U.
The optimization problems that involve the mean squared error loss function can be solved using the simplex method for quadratic programming. The solution in this case is unique. The uniqueness can be observed geometrically, i.e., the optimal solution is the point of contact of the convex simplex formed by the linear constraints and the nonlinear strictly convex cone formed by the objective function. Convex and strictly convex shapes can touch each other in only one point.
In the appendix of this paper, we investigate if the linear optimization problems can be solved using non-analytical (combinatorial) approaches. Namely, the dual versions of problems (12) and (13) can be modeled as the min-cost flow problem and its variations. We recall the min-cost flow problem and some algorithms used to solve it in Appendix A while we show how to model dual problems of (12) and (13) as the min-cost flow problem and a variation of the min-cost flow problem, respectively, in Appendix B. In the same section, we provide a greedy algorithm to solve the aforementioned variation based on the algorithm that solves the original min-cost flow problem. Since the algorithm is new, we provide its proof of correctness in Appendix C.

Properties
In this section, we prove some properties of the granular approximations obtained in Section 4. First we show that the proposed approach is a generalization of the KS approach when the number of ordinal labels is 2.
Proof. IfĀ is crisp, it is obvious that the objective function from (11) corresponds to the objective function from (5) for K = 2, where the labels with value 1 are those that are more preferred. Regarding the constraints, we examine the granular conditions in form R(u, v) ≤ I(Â(v),Â(u)). If R(u, v) = 0, then there are no restrictions on the implication, i.e., we do not have a constraint. If R(u, v) = 1 thenÂ(v) ≤Â(u) from the ordering property of I (8). Since R(u, v) = 1 means that u v (u dominates v) then the condition R(u, v) = 1 ⇒Â(u) ≥Â(v) is equivalent to u v ⇒Â(u) ≥Â(v) which is exactly the condition from (5).
In the remainder of this section, we prove some properties of our approach that also hold in the original (crisp) one.
For granular approximations obtained with the p-quantile loss, the monotonicity property holds.
Proposition 5.2. Let p and q be two real numbers from the unit interval and letÂ p andÂ q be the outputs of the optimization problem (12) or (13) with p and q as probability parameters. It holds that Proof. The proof is provided in Appendix D. It relies on the greedy combinatorial approach presented in the previous sections of Appendix, hence those previous sections are necessary for the understanding of the proof.
Before we introduce the next property, we present the following lemma. Proof. For A being granularly representable, we have Applying negation N on both sides of the inequality, we have The first equivalence follows from Proposition (2.4) while the second is the residuation property.
In the previous proposition, implication becomes an equivalence if we use IMTL triplets as operators.
For the fuzzy rough approximations that are obtained with IMTL operators, we have the well known duality property as stated in Subsection 2.3. That property can be extended to the granular approximations. Proposition 5.3. Let α * u , u ∈ U be a minimizer of the optimization problem (11) with nilpotent t-norm T L,ϕ , relation R, observationsĀ and risk u∈U L p (ϕ(Ā(u)), α u ) (for short L p problem). Then 1 − α * u , u ∈ U, is a minimizer of the optimization problem (11) with the same t-norm, relation R −1 , observationsĀ and risk u∈U L 1−p (ϕ(coĀ(u)), α u ) (for short L 1−p problem). Proof. Let α u , u ∈ U, be a feasible solution of the L p problem, i.e., it satisfies granularity conditions w.r.t. relation R The expression above is equivalent to which states that 1 − α u , u ∈ U satisfies granularity conditions w.r.t. relation R −1 and, therefore, it is a feasible solution of the L 1−p problem. We observe that ϕ(coĀ(u)) = ϕ(ϕ −1 (1 − ϕ(Ā * (u)))) = 1 − ϕ(Ā * (u)). Regarding the loss function, we have that Since α u , u ∈ U, is a feasible solution of the L p problem if and only if 1 − α u , u ∈ U, is a feasible solution of the L 1−p problem, and due the previous equality, we finish the proof.
The duality also holds for the mean squared error risk. The proof is very similar to the proof of Proposition 5.3 where the only difference is that the loss function stays the same in the dual problems.
The next step is to examine fuzzy rough lower and upper approximations in the context of the newly defined optimization problems. In the crisp case, the lower and upper fuzzy rough approximations are seen as sets of necessary and possible knowledge respectively. In other words, the actual ill-known knowledge must contain the lower approximation and be contained in the upper one. In probabilistic terms, the probability that the actual knowledge is between these approximations is 1 [47]. Hence, the lower and upper approximations are the extreme values in the probability distributions of the actual knowledge. It means that the lower approximation is the 0-quantile while the upper approximation is the 1-quantile. We have the following proposition. Proof. When optimization problems (12) and (13) are considered in terms of α and not in terms ofÂ, they can be seen as problem (11) with t-norm T L or T P , relation R ϕ and observationsĀ ϕ . If p = 1, then the loss function for u ∈ U is equal 0 if α u −Ā ϕ (u) ≥ 0 and to a positive value otherwise. If for all u ∈ U it holds that α u ≥Ā ϕ (u), then the objective is 0 and hence any such α is a solution. Such fuzzy set α contains fuzzy setĀ ϕ and is granularly representable w.r.t. t-norm T L or T P and relation R ϕ . From Proposition 2.5, the smallest such α is the fuzzy rough upper approximation, i.e., the smallest solution is or with T P instead of T L . Then, the final solutionÂ * is obtained The derivation for T P is the same.
The proof for the lower approximation is analogous.

Didactic Examples
In this section, we provide two didactic examples of the approaches (12) and (14) that use T L,ϕ t-norms. One example is related to a binary classification problem while the other one is related to a regression problem.
To evaluate pairwise fuzzy relation values, in both cases we use a T Lequivalence relation called triangular similarity. For a condition attribute q, it is defined as while the overall relation is then R(u, v) = min q R q (u, v). More details on such similarity relation are provided in [38]. For the classification purpose, we use 4 instances from two classes of the well-known iris dataset which can be found in the UCI dataset repository [48]. Those instances are shown in Table 2.  In this binary classification example, the observed values are from the set {0, 1} where they indicate if an instance is in a decision class or not. We consider such crisp set as a fuzzy one in order to obtain granular approximations. After applying the relation on every pair of instances, we get a relation matrix shown in Table 3.  In Table 3, the first row and the first column stand for indices of instances from Table 2. The remaining entries are values of the relation on the corresponding pair of instances.
After the matrix is calculated, we pass it together with the decision attribute to the optimization problem (12)    In every row, we have a granular approximation for a corresponding probability parameter from the first column. Every entry is a fuzzy membership degree for the corresponding instance which may be interpreted as the degree up to which the instance belongs to class with label 1. Since that fuzzy value is unknown, we have its distribution characterized with quantiles. For example, in the second row of Table 4, we say that with probability 0.25, the degree up to which instance 3 belongs to the class with label 1 is not greater than 0.475.
The results for the mean squared error used in (14) are shown in Table 5.  In this case, we may say that the expected degree to which instance 3 belongs to the class with label 1 is 0.696.
In this particular case, the granularity (consistency) condition expressed in form T ( R(u, v),Â(v)) ≤Â(u) can be interpreted as "if instance v belongs to the decision class, and u is similar to v, then u also belongs to the same decision class". We express it numerically for instances 1 and 3, for which the similarity can be obtained from Table 3 as 0.525. We notice that the granularity condition was not satisfied with their original labels from Table 2 since it holds that T (0.525, 1) = 0.525 > 0. Using the new labels from Table  5, we have that the condition is satisfied since T (0.525, 0.696) = 0.221 ≤ 0.221 and T (0.525, 0.221) = 0 ≤ 0.696.
In the regression example, we use 'Real estate valuation' dataset that can also be found in the UCI dataset repository [48]. The goal of this dataset is to predict a price of a real estate given its features like position, age, closeness to some important facilities, etc. We use 5 instances displayed in Table 6. Since the decision values are real estate prices, we have to model them as a fuzzy set in order to use our methods. Fuzzy predicate "expensive" can be modeled using prices, hence we construct a new fuzzy set which indicates how expensive is a certain real estate based on the prices. The new fuzzy set is constructed using the following formula: where Q stands for p-quantile. This is similar to the linear scaling to the [0, 1] range, just here using quantiles we achieve that 0.5 percent of smallest prices has expensiveness 0 and 0.5 percent of largest prices has expensiveness 1. We use this approach in order to handle extreme values -outliers. After applying this transformation on the complete dataset, the fuzzy values of our chosen 5 instances are given in the last column of  Applying the relation on every pair of instances leads to the relation matrix shown in Table 7.  As before, we pass matrix values with the new expensiveness values to optimization problem (12) with probability parameters p ∈ {0, 0.25, 0.5, 0.75, 1}. The obtained granular approximations are given in Table 8.  The obtained fuzzy values are estimations of quantiles of the expensiveness, under the assumption that it is a random fuzzy set and that its realiza-tions are given in Table 6. We interpret the values in a way that, for example, in the third row of Table 8 we say that the expensiveness of instance 2 is less than 0.54 with probability 0.5, or in the fourth row of the table, we say that the expensiveness of instance 4 is less than 0.843 with probability 0.75.
The results for the mean squared error used in (14) are shown in Table 9.  In this case, we may say that the expected expensiveness of instance 4 is equal to 0.695.
In the regression case, the granularity (consistency) condition expressed in form T ( R(u, v),Â(v)) ≤Â(u) is interpreted as "if instance v is expensive and u is similar to v, then u is also expensive". We observe instances 2 and 4 for which the similarity is equal to 0.696. Their original labels from Table 6 do not satisfy the granularity condition since T (0.696, 0.938) = 0.634 > 0.54. Using the new labels from Table 9, we have that T (0.696, 0.695) = 0.391 ≤ 0.54 and T (0.696, 0.54) = 0.236 ≤ 0.695. Hence, the condition is satisfied.

Conclusion and future work
In this paper, we introduced a novel machine learning approach for handling inconsistencies in prediction problems with respect to a fuzzy relation. Our work was motivated by the method introduced by Kot lowski and S lowiński [10] for handling monotone inconsistency and we showed that the novel approach is a generalization of the same method in the binary classification case. Using fuzzy relations, the novel method is able to handle gradual relationships among instances while the KS approach can distinguish only two cases: either instances relate or not.
The novel approach produces a granular approximation of a fuzzy set. The approximation is granularly representable (without inconsistencies) and as close as possible to the original fuzzy set (w.r.t. a given loss function). It can be seen as a fuzzy counterpart of the monotone approximation produced by the KS approach. As in the work of Kot lowski and S lowiński, we provided statistical foundations of the granular approximations. In the next step, we formulated optimization problems in order to calculate the approximations and we showed some of their important properties. At the end, we provided two didactic examples; one for a binary classification problem and one for a regression problem. In the didactic examples, we showed how fuzzy relations are used to model relationship among numerical data, how the granular approximations are calculated and how to interpret them in the two cases for different loss functions.
Since this work is theoretical in its nature, for the future work we propose its implementation and verification on practical applications. As it was already mentioned in the Introduction, the possible applications are in fuzzy rough set based methods and fuzzy rule induction. Therefore, our future work will mainly concentrate on the aforementioned implementation and computational experiments.

Appendix A. Minimum-cost flow problem
This section is based on the monograph [49], especially on its 9th chapter. A flow network is defined as a directed graph where a real value called imbalance is assigned to each node. Imbalances split nodes into two subsets: supply nodes with a positive imbalance (supply value) and demand nodes with a negative imbalance (demand value). Moreover, each edge is characterized by a positive real capacity, and a cost value. We also assign flow amounts to each edge which satisfy the condition that they are at most as large as capacities. More formally, let G be a finite set of nodes, E ⊆ G × G the finite set of edges, while F = (G, E) is the flow network. We denote imbalances with b i for i ∈ G, capacities with l i,j , costs with c i,j and flow with z i,j for (i, j) ∈ E.
The minimum-cost flow problem is an optimization problem defined on a flow network where we want to transport flow from the supply nodes to the demand nodes, such that -the difference between the flow that leaves a node and the flow that enters the node is equal to the imbalance of this node, -a flow in a particular edge is at most as large as the capacity of that edge, and -the total cost of the flow transportation is minimal.
Formally, we have the following problem: We distinguish two sets of constraints in the previous optimization problem: balance constraints (A.1b) and capacity constraints (A.1c). If we sum the balance constraints, we get i∈G b i = 0 which states that the amount of supply is equal to the amount of demand, which is a necessary assumption to have a feasible solution.
We say that a flow is feasible if it is a feasible solution of (A.1), while we say that we have a pseudo-flow if only the capacity constraints are satisfied.
For a given pseudo-flow z ′ , a residual network F ′ = (G, E ′ ) can be defined. We have new imbalances: while for each edge (i, j) ∈ E for which z ′ i,j > 0, we add the reverse edge (j, i) to the network with cost c ′ j,i = −c i,j , while keeping the original edge. The capacity of the original edge (i, j) in F ′ is l ′ i,j = l i,j −z i,j , while the capacity of the added reverse edge (j, i) is l ′ j,i = z i,j . We may notice that when adding a new edge (j, i) to E ′ , there can already exist an edge (j, i) from E. However, in our case of use, we will not face such an issue, i.e., we will have either (i, j) or (j, i) in E and not both at the same time. The residual network keeps the complete information about flow z ′ which can be reconstructed from F ′ .
The concept of residual network is important for the development of algorithms for solving (A.1). In this moment, we will not discuss the existence of a feasible solution in general since later we will show that it always exists in our case of use.
A cost of a particular path or cycle in the flow network is calculated as the sum of the costs of edges in that path or cycle. For an optimal flow z * , we have the following result.
Proposition Appendix A.1. A flow z * is optimal if and only if there are no cycles of negative cost in the residual network F (z * ).
Bearing in mind Proposition Appendix A.1, a simple algorithm can be constructed to solve (A.1). Namely, we construct an initial feasible flow in our network, then search for the negative cycles and eliminate them.
However, a more useful algorithm for us is the Successive Shortest Path (SSP) algorithm for solving the minimum-cost flow problem. The algorithm is provided as Algorithm 1. Pick supply node i and demand node j 7: Calculate the shortest path P from i to j using cost values from F ′

8:
Send the largest possible amount of flow through P 9: Update F ′ 10: Reconstruct z from F ′ The shortest path P can be calculated using the Bellman-Ford algorithm since F ′ may contain negative values. The largest possible amount of flow through P is calculated as The residual network is then updated such that ∈ P The idea of the proof of correctness is that sending a flow through the shortest path does not produce negative cycles in the residual network. Hence, when all supply is sent to the demand nodes and the feasible solution is achieved, it will be an optimal one.
We also introduce generalized network flows based on Chapter 15 of [49]. In some cases, the flow in a particular edge may be increased or decreased by a multiplier after it leaves the left node of the edge. Denote the multipliers with m i,j for (i, j) ∈ E. The generalized minimum-cost flow problem is then formulated as If the multiplier is greater than 1, then the flow is increased while if it is smaller than 1, then the flow is decreased.
Different theoretical results hold for the generalized minimum-cost flow problem (A.2). Fortunately, our particular case of (A.2) allows obtaining similar properties as we have in the ordinary minimum-cost flow problem (A.1).

Appendix B. Duality and the combinatorial approach
In this section, the dual optimization problems of (12) and (13) are considered. In our particular case, the dual problems are interesting since they can be modeled using graph theory and can be solved using combinatorial optimization methods. These combinatorial algorithms may not be more efficient than the simplex method used for solving linear programs, but their development is important since they allow us to prove some interesting properties of the estimated fuzzy set. We examine optimization problem (12). First, we eliminate variables x u , u ∈ U, using constraints x u = y u +Ā ϕ (u) − α u and we denote M(u, v) = 1 − R ϕ (u, v). Then, the problem is reformulated as (B.1) Its dual problem is then In (B.2), variables z u,v , u, v ∈ U, correspond to the first set of constraints from primal (B.1), while variables z 0,u , u ∈ U, correspond to the second set of constraints from the primal. The first set of constraints in (B.2) corresponds to variables α u , u ∈ U, from the primal, while the second set of constraints corresponds to variables y u , u ∈ U, from the primal.
If we sum up the equality constraints, we get u∈U z 0,u = np where n = |U|. Bearing that in mind, we see that (B.2) is exactly the minimum-cost flow problem on n + 1 nodes where we have one supply node with imbalance b 0 = np and n demand nodes with imbalances −p. From the supply node, to all other nodes we have flow z 0,u , costsĀ ϕ (u), while all capacities are equal to 1. Among the demand nodes, there is a flow z u,v , u, v ∈ U, costs M(u, v), and there are no capacity constraints.
To make our model even simpler, we utilize the T -transitivity of the relation R. It is easy to verify that the T -transitivity is equivalent to M(u, v) + M(v, w) ≥ M(u, w) for u, v, w ∈ U. Using this fact, we have that there is an optimal flow which does not use two consecutive edges that are between demand nodes. Assume that for an optimal flow z * we have z * u,v > 0 and z * v,w > 0, and let δ = min(z * u,v , z * v,w ). Then the flow z * u,v − δ, z * v,w − δ, z * u,w + δ is feasible and at most as expensive as the previous flow, i.e., it is optimal. The new flow does not use two consecutive edges since either z * u,v − δ or z * v,w − δ is 0. The previous elaboration further implies that an optimal flow from the supply node can travel through at most one intermediary node to the destination demand node. Hence, our initial network flow on n + 1 nodes can be transformed into a flow network on 2n + 1 nodes which has a form of a bipartite graph plus the supply node. One independent set in the bipartite graph is formed by the intermediate nodes, while the other independent set is formed by the destination nodes. The bipartite network representation is useful from the perspective of the flow decomposition. For a feasible flow, it is easy to represent it as a sum of simple flows that go from the supply node to the destination node. In the original network, one node can be a destination node for some flow but also an intermediate node for a different flow. Hence, the decomposition is harder in the original network. The decomposition will be important later when dealing with the dual of (13).
The next question is how to reconstruct optimal solution of the primal problem, i.e., to calculate α * from a solution of the dual z * . Following the duality theory provided in [49], an optimal vector α * can be obtained as lengths of shortest paths from the supply node to the corresponding destination nodes in the residual network of z * . Now, we examine the dual of (13). The linear program here can be rewritten similarly as (B.1), just with the different granularity constraints.
The difference between (B.2) and (B.3) is that in the latter, we have multipliers R ϕ (u, v), u, v ∈ U, instead of costs on the edges. Due to the multipliers, we now deal with the minimum-cost flow problem on a generalized flow network with n + 1 nodes among which there are n demand nodes with demand −p and one supply node with an unspecified amount of supply.
We may notice that in this case the edges of the network consist of two different groups. The first group is formed by the edges from the supply nodes to the demand nodes. Those edges have costs and do not have multipliers. The second group is formed by the edges among the demand nodes. Those edges, conversely, have multipliers and do not have costs. Similarly to (B.2), we are able to utilize the T -transitivity of R ϕ w.r.t. T P in a way that there is an optimal flow which does not use two consecutive edges from the second group. If we have three demand nodes u, v, w ∈ U in a network and an optimal flow that uses edges (u, v) and (v, w), we can redirect the flow to use only edge (u, w) and the redirected flow will have smaller or equal loss than the original flow. That will further lead to the smaller or equal cost of the redirected flow which makes it optimal. Therefore, as above, there is an optimal solution in which a flow travels from the supply node to the destination demand node using at most one intermediate node. That again further implies that the initial general network on n + 1 nodes can be transformed into a generalized bipartite flow network on 2n + 1 nodes. For the new network, the same model applies as in Figure B.1. Using this model, we can clearly see the difference between two groups of edges introduced above. The first group is formed by the edges between the supply node and the left partition of the bipartite graph (intermediate nodes), while the second group is formed by the edges between the two partitions of the bipartite graph.
As before, for a given flow on the generalized bipartite network, we have the corresponding residual network. The same properties apply as above except the case when the flow passes through an edge with multiplier. In that case, if the original edge has multiplier R ϕ (u, v) then the reverse edge in the residual network will have multiplier 1 Rϕ(u,v) which is an edge of a gain type (greater than 1).
We will now construct a new algorithm for solving a generalized minimumcost flow problem on a generalized bipartite flow network. The algorithm is based on the existing SSP algorithm presented in Algorithm 1. Assume that we have a demand node f u to which we want to deliver some flow b. We want to deliver the flow at the cheapest possible price. If we deliver a flow using intermediate node e v , then the amount of flow that we have to take from the supply node is b Rϕ(v,u) and the cost of such flow is bĀϕ(v) Rϕ(v,u) . In general, a price to deliver a unit of flow is a ratio of the cost of an edge from the supply to the first partition and the product of multipliers of edges that connect the two partitions. Bear in mind that in the residual network, a flow may use multiple edges between partitions (edges with multipliers) to deliver the flow. Using this, we construct the greedy approach presented as Algorithm 2. Calculate the smallest possible cost from the supply node to i 8: Calculate the largest amount of flow that can be sent through the least costly path 9: Send the calculated flow through the least costly path 10: Update F ′ 11: Reconstruct z ′ from F ′ To calculate the smallest possible cost from the supply node, we can use a shortest path method. We want to minimize the ratio of one cost value (from the supply to the intermediate nodes) and a product of multipliers (between intermediate and destination nodes). If we apply logarithms on the cost values and reciprocals of the multipliers, we may apply the Bellman-Ford algorithm to calculate the shortest path between the supply node and the chosen demand node in order to obtain a least costly way to transport the flow.
After the shortest path is determined, we have to calculate the amount of flow that will be taken from the supply node in order to deliver the maximal amount of flow to the demand node. In comparison with the standard minimum-cost flow problem, here we have to take into account all the loses and gains that happen during the flow transfer. Denote the shortest path in the residual network with P = (0, e u 1 , f u 2 , e u 3 , . . . , f u k ) and let b be a demand of node f u k . We would like to deliver |b| (| · | stands for absolute value) amount of flow to the demand node from the supply node, but that is not always possible due to the capacities of particular edges on path P . The maximal amount of flow can be determined recursively. The maximal amount of flow that can be transferred from node f u k−2 to node f u k is bounded with the capacity of the reverse edge l ′ fu k−2 ,eu k−1 and the demand divided with the rational representation a = q r for q and r being integers. Let C be the least common multiple (LCM) of all integers q and r for all parameters in (B.3). If we multiply the right side of the constraints in (B.3) with C, then all the demand value will become integers and all intermediate flows in Algorithm 2 will become integers. That further implies that all the updates on demands in Algorithm 2 will be integer which further implies that the algorithm will terminate in at most Cpn steps.
In practice, the termination is always guaranteed since computers can work only with rational numbers Now, lets define a flow cycle in the residual generalized bipartite network. The cycle starts with an edge from the first part (costly edges without multipliers) of the network, then it contains edges from the second part (edges with multipliers without costs) and ends with a reverse edge from the first part. A model of such cycle is shown in Figure C In Figure C.2, the dashed line between e u 1 and e u 2 stands for the subpath that contains only the edges from the second part of the residual network. Also, it may hold that e u 1 ≡ e u 2 . In that case, the cycle consists only of the edges from the second part. Let R ϕ (e u 1 , e u 2 ) be a multiplier of the path that consists of the edges from the second part of the residual network, i.e., a product of the multipliers on the edges from the path. We say that the cycle is of negative cost if A ϕ (u 1 ) < R ϕ (e u 1 , e u 2 )A ϕ (u 2 ). For the reminder, A ϕ (u 1 ) and A ϕ (u 2 ) are the costs on edges (0, e u 1 ) and (0, e u 2 ). The reason why the cycle is of negative cost is that if we send a unit of flow along it, the cost of that flow is A ϕ (u 1 ) − R ϕ (e u 1 , e u 2 )A ϕ (u 2 ), i.e., the cost is negative. Such flow would not change any demand value on the destination nodes but it will reduce the overall cost of the flow.
The next proposition utilizes the bipartite representation of the flow network.
Proposition Appendix C.2. Every flow in a generalized bipartite network can be represented as a sum of a finite number of simple path flows from the supply node to a destination node.
Proof. Let z be a flow and consider an edge (e u 1 , f u 2 ) from the second part of the network and its flow value z eu 1 ,fu 2 . That edge receives a flow from edge (0, e u 1 ) which is a part of path flow z P from path P = (0, e u 1 , f u 2 ) that connects the supply node and the destination node f u 2 . z P is then a summand in the representation while the remaining flow z − z P has no flow on the edge (e u 1 , f u 2 ) and hence we can remove that edge from the network flow. If we continue, in every step we will construct one summand and remove one edge from the second part of the network. Since we have a finite number of edges, we have a finite number of summands.
We have the following result.
Proposition Appendix C.3. Solution z * is optimal in the generalized bipartite network if and only if its residual network does not contain negative cost cycles.
Proof. (⇒) When the solution is optimal, there are no negative cost cycles. If otherwise, we could send a flow through a negative cost cycle and we would decrease the cost of the overall flow as described above. That contradicts the optimality.
(⇐) Assume that z * is a feasible solution whose residual network does not contain negative cost cycles and let z ′ be a feasible solution. Let z ′ = z * + z ′′ . We first show that z ′′ is a feasible flow from the residual network of z * , i.e., it satisfies its constraints. For an edge (0, e u 1 ) if the flows are different, we can have either z ′ 0,eu 1 > z * 0,eu 1 or z ′ 0,eu 1 < z * 0,eu 1 . In the first case, it holds that z ′ 0,eu 1 = z * 0,eu 1 + z ′′ 0,eu 1 , i.e., z ′′ 0,eu 1 uses the original edge. Since z ′ 0,eu 1 ≤ 1 then z ′′ 0,eu 1 ≤ 1 − z * 0,eu 1 which is a constraint from the residual network. In the second case, it holds that z ′ 0,eu 1 = z * 0,eu 1 − z ′′ eu 1 ,0 , i.e., z ′′ eu 1 ,0 uses the reverse edge. Since z ′ 0,eu 1 ≥ 0 then z ′′ eu 1 ,0 ≤ z * 0,eu 1 which is a constraint for the reverse edge from the residual network. Using similar reasoning, we can conclude the same for the whole network.
The next step is to show that z ′′ is a sum of a finite number of simple flow cycles, as shown in Figure C.2, i.e., it has a cycle representation. Proposition Appendix C.2 states that both flows z ′ and z * are sums of simple flows on paths from the supply node to a destination node. Take a summand z ′ P 1 of z ′ and summand z * P 2 of z * for P 1 = (0, e u 1 , f u 3 ) and P 2 = (0, e u 2 , f u 3 ). The paths have the same destination node. Assume that the first summand delivers amount b 1 of flow to the destination node while the second delivers amount b 2 of flow to the same node. W.L.O.G. assume that b 1 ≥ b 2 . Then the flow b 2 b 1 z ′ P 1 − z * P 2 is a flow along cycle (0, e u 1 , f u 3 , e u 2 , 0) and one of the summands in the cycle representation of z ′′ . After the summand is identified, we remove its flow from the consideration. In that moment, z * P 2 is fully removed while we are left with (1 − b 2 b 1 )z ′ P 1 from the first path. We continue to create flow cycles as summands from the remaining path flows from z ′ and z * . Since after every summand is identified we remove one path flow, the number of summands is finite. Hence, z ′′ is a sum of a finite number of cycle flows. Since z ′′ is a flow in the residual network of z * , all the cycles from its cycle representation are of positive cost by the assumption which implies that z ′′ is of positive cost. Since the cost of z ′ is a sum of costs of z * and z ′′ , cost of z ′ is larger than the cost of z * . Since flow z ′ was an arbitrary feasible flow, we conclude that z * is an optimal flow. Proposition Appendix C.4. Algorithm 2 returns an optimal solution.
Proof. Assume that in one iteration of Algorithm (2), the shortest path had the form P 1 = (0, e u 2 , . . . , f u 3 ) and that after the step, the negative cost cycle (0, e u 1 , . . . , f u 3 , . . . , e u 2 , 0) was formed. The negative cost cycle is formed from the path P 2 = (0, e u 1 , . . . , f u 3 ) and the reverse path P 1 . The model of such cycle is represented in Figure C The dots in the cycle as well as dashed edges in the figure stand for edges from the second part of the residual network (edges with multipliers). If the cycle is negative, then A ϕ (u 1 ) < R ϕ (e u 1 , f u 3 ) R ϕ (f u 3 , e u 2 )A ϕ (u 2 ). The latter is equivalent to Aϕ(u 1 ) Rϕ(eu 1 ,fu 3 ) < Aϕ(u 2 ) Rϕ(eu 2 ,fu 3 ) which states that path P 2 is actually shorter than P 1 which contradicts the assumption that P 1 is the shortest path at this step. Hence, at every iteration of Algorithm 2, there are no negative cost cycles and as soon as the feasible solution is achieved, it will be an optimal one according to Proposition Appendix C.3.
After we constructed the algorithm that solves the dual optimization problem, we need to obtain an optimal solution for the primal which was our initial goal. First, we need one technical proposition.
Proposition Appendix C.5. For a given generalized bipartite network, there exists optimal solution z * for which it holds z * 0,eu > 0 =⇒ z * eu,fu > 0. Proof. Assume that for some solution z * and some instance u we have that z * 0,eu > 0 and z * eu,fu = 0. Then, in the simple path decomposition of the flow, we have path (0, e v , f u ) that delivers flow to f u , and path (0, e u , f w ) that uses flow from edge (0, e u ). Then, in the residual network of z * , C = (e u , f u , e v , f w , e u ) is a cycle. Due to transitivity of R, it holds R ϕ (v, u) R ϕ (u, w) ≤ R ϕ (v, w).
If R ϕ (v, u) R ϕ (u, w) < R ϕ (v, w), then C is a negative cost cycle which contradicts the optimality of z * . If R ϕ (v, u)R ϕ (u, w) = R ϕ (v, w) then cycle C is a zero-cost cycle and a flow can be sent through the cycle without violating optimlaity. Hence, sending some amount of flow through the cycle, we will construct another optimal solution z * * where z * * eu,fu > 0. In practice, if we obtain an optimal solution containing an edge for which the previous proposition does not hold, we can get another optimal solution, without such edges, as explained in the proof of the previous proposition. From now on, we assume that we have an optimal solution for which the previous proposition holds.
We continue with the duality theory of the linear programs. According to the strong duality theorem [50], if there exists an optimal solution of the dual problem z * then, there exists an optimal solution for the primal problem α * , and it holds that the values of objectives in (B.1) and in (B.3) are equal, i.e., In the previous expression, y u is replaced with its definition. In an optimal solution, we have that • z * u,v > 0 =⇒ α * u − R ϕ (u, v)α * v = 0, for u, v ∈ U. We have the following conclusion: if we solve the dual optimization problem and obtain an optimal solution z * , then a solution of the primal optimization problem is any α * which satisfies the conditions listed above.
Moreover, α * can be constructed by performing step 7 of Algorithm 2 on the residual network of z * , i.e., it is the smallest possible cost of the transport from the supply node to the destination nodes. It is easily verifiable that such α * satisfies the conditions above. The proof of this verification lies in that if we assume that some condition is not satisfied, then we would have a negative cost cycle which contradicts the optimality of z * . To prove the contradiction, we need Proposition Appendix C.5.
the original edges between the same nodes will not reduce the costs since the shortest paths now chose among the smaller set of edges. The same holds if we remove reverse edges between the intermediate nodes.
The last step is to prove that adding reverse edges between the intermediate and destination nodes will not reduce the costs from the supply to the destination nodes.
For that purpose, we consider Figure D With dashed lines, we denote certain paths whose costs are marked on the figure. In both cases of T L and T P , the costs are the values used to calculate the shortest paths. Assume that in step i, we were calculating the shortest path between 0 and f u 2 and we obtained that the shortest path is (0, . . . , e u 1 , f u 1 , . . . , f u 2 ) and since some flow is sent through that path, a reverse edge (f u 1 , e u 1 ) is created with cost −x. Assume that before step i, the shortest path from 0 to f u 3 was (0, . . . , f u 3 ) with cost C y while after the previous step and after adding reverse edge (f u 1 , e u 1 ) the shortest path is (0, . . . , f u 1 , e u 1 , . . . , f u 3 ) with cost C a − x + C b . Then, we have that C a + C b < x + C y . Since the shortest path in step i was (0, . . . , e u 1 , f u 1 , . . . , f u 2 ), it holds that C c + x ≤ C a . Adding this to the previous expression, we have that The last inequality contradicts the assumption that before step i, the smallest cost between 0 and f u 3 is C y .