Centered solutions for uncertain linear equations

Our contribution is twofold. Firstly, for a system of uncertain linear equations where the uncertainties are column-wise and reside in general convex sets, we derive convex representations for united and tolerable solution sets. Secondly, to obtain centered solutions for uncertain linear equations, we develop a new method based on adjustable robust optimization (ARO) techniques to compute the maximum size inscribed convex body (MCB) of the set of the solutions. In general, the obtained MCB is an inner approximation of the solution set, and its center is a potential solution to the system. We use recent results from ARO to characterize for which convex bodies the obtained MCB is optimal. We compare our method both theoretically and numerically with an existing method that minimizes the worst-case violation. Applications to the input–output model, Colley’s Matrix Rankings and Article Influence Scores demonstrate the advantages of the new method.


Introduction
Systems of linear equations are of immense importance in mathematics and its applications in physics, economics, engineering, and many more fields. However, the presence of unavoidable errors (inaccuracies) in the specification of parameters in both the right- and left-hand sides introduces uncertainty in the sought solution. The uncertainties may be raised due to measurement/rounding errors in the data of the physical problems, estimation errors in the estimated parameters by using expert opinions and/or historical data, or numerical errors associated with finite representation of numbers by computer (see Ben-Tal et al. 2009).
A basic version of the problem that we consider is well known in the context of interval linear systems. For a given system of linear equations, where the coefficient matrix A ∈ R m×n and right-hand side b ∈ R m are uncertain and allowed to vary uniformly and independently of each other in the given intervals: where a i j , a i j , b i , b i ∈ R, for all i, j, are the lower-and upper-bounds of the components in the matrix A and vector b, respectively. Let U = A×B. There are three well known solution sets to the linear system (1), see, e.g., , Shary (2011), Popova (2012) and Hladík and Popova (2015): united solution set controllable solution set tolerable solution set For a united solution set X , each of the possible (A, b) ∈ U has equal claim to be the true realization of the physical problem. Different parameters may produce different solutions for the system of linear equations. Since the pioneer work by Oettli and Prager (1964), much literature has been devoted to describe the ranges of the components of the solution x ∈ X for interval linear systems, i.e., where x i denotes the i-th element of vector x. The main source of difficulties connected with obtaining ranges of x i is the complicated structure of the solution set X . Although the intersection of the solution set and each orthant is a convex polyhedron, the union of those polyhedra, i.e., X , is generally non-convex. Oettli (1965) proposes to use a linear programming procedure in each orthant (i.e., 2 n orthants in total) to determine x i and x i , for all i ∈ [n]. Rohn and Kreinovich (1995) show that, in general, determining the exact ranges for the components of x ∈ X is an NP-hard problem. For a comprehensive treatment and for references to the literature on interval linear systems one may refer to the books Neumaier (1990), , Fiedler et al. (2006) and Moore et al. (2009). Due to the NP-hardness of solving (3) exactly, many ingenious methods have been developed to obtain sufficiently close outer estimates of the solution set, e.g., Hansen (1992), Jansson (1997), Ning and Kearfott (1997), Rump (2010), Hladík (2014), Popova (2004Popova ( , 2014, Rohn (1981), Alefeld et al. (1998), Calafiore and Ghaoui (2004), Popova and Krämer (2007) and Hladík (2012) consider interval linear systems with dependent data. We refrain here from listing papers dedicated to computing enclosures since they are simply too many. Interval linear systems has been applied to many engineering problems described by systems of linear equations involving uncertainties. These problems include analysis of mechanical structures (see Smith et al. 2012;Muhanna and Erdolen 2006), electrical circuit designs (see Dreyer 2005;Kolev 1993) and chemical engineering (see Gau and Stadtherr 2002). For more applications we refer to the book Moore et al. (2009).
In this paper, we consider the system of uncertain linear equations: where the coefficient matrix A : R n ζ → R m×n and right-hand side b : R n ζ → R m are affine in ζ , and the uncertain parameters ζ ∈ R n ζ reside in the uncertainty set U. Firstly, we focus on systems of linear equations with column-wise uncertainties. Let A(ζ ) = [a ·1 (ζ 1 ) a ·2 (ζ 2 ) · · · a ·n (ζ n )], b(ζ ) = b(ζ 0 ), and ζ = [ζ 0 ζ 1 · · · ζ n ] . We represent the system (4) as follows: where the vector function a · j is affine in ζ j ∈ U j , the vector function b is affine in ζ 0 ∈ U 0 , and the set U j is convex, for all j ∈ [n] ∪ {0}. The corresponding united solution set X is: where U = × n j=0 U j , and ζ = [ζ 0 ζ 1 · · · ζ n ] ∈ R n ζ . Most of the existing literature consider (4) with independent interval uncertainties. Oettli (1965) shows that for a special class of (6), where uncertainties are component-wise and reside in independent intervals, i.e., (2), the intersection of the united solution set X and any orthant of R n is a convex polyhedron. Popova (2014) devises a method to obtain close outer estimates of (4) with independent interval uncertainties. Popova (2015) characterize the solvability of (4) with independent interval uncertainties. Here, we consider uncertain linear systems with column-wise (dependent) uncertainties that reside in general convex sets, and derive convex representations of X in an arbitrary orthant. Moreover, when U j are polyhedral, j ∈ [n] ∪ {0}, the set X is also polyhedral; when U j are ellipsoidal, j ∈ [n] ∪ {0}, the set X is conic quadratic representable.
Moreover, via a convex representation technique and the techniques from (adjustable) robust optimization, we are able to derive convex representations of a special class of X con . We further show that X tol can be interpreted as a solution set for a static robust optimization problem, which admits tractable reformulations for many important classes of A and B.
We propose to compute the maximum size inscribed convex body (MCB) of the set of possible solutions for systems of uncertain linear equations. It is intuitively appealing to find a centered solution that is "far" from the boundaries of the solution set X (i.e., infeasibility). The obtained MCB is an inner approximation of the solution set, and its center is a potential solution to the system. We extend the method developed in Zhen and den Hertog (2017) to compute the MCB of the solution set X . This method is independently developed in Zhang et al. (2017) and applied to optimal control problems. Furthermore, we use the recent results of Zhen et al. (2016) for adjustable robust optimization (ARO) to characterize for which convex bodies the obtained MCB is optimal.
We compare our new method both theoretically and numerically with an existing method. The conventional method of determining a robust solution for systems of uncertain linear equations first appears in the context of robust least-squares (RLS) problems Ghaoui and Lebret (1997). The RLS method finds the minimizer x RL S of the worst-case 2-norm violation of the system: The tractability of Problem (7) strongly relies on the choice of the uncertainty set U. Ben-Tal et al. (2009) show that Problem (7) under independent interval uncertainties can be reformulated into an SOCP problem. In Ghaoui and Lebret (1997), Beck and Eldar (2006) and Jeyakumar and Li (2014), authors derive an SOCP or a semidefinite programming (SDP) reformulation of Problem (7) under ellipsoidal uncertainties. Burer (2012) and Juditsky and Polyak (2012) solve (7) to find the robust rating vectors for Colley's Matrix Ranking and Google's PageRank, respectively. The contributions of this paper may be summarized as follows: 1. For column-wise (dependent) uncertainties that reside in convex sets, we show that the united solution set X in any orthant is convex, and derive a convex representation of X . We then derive a convex representation of X con via robust optimization techniques, and further discuss the cases in which X tol can be reformulated into equivalent convex sets. 2. Based on ARO techniques, we develop a new method to compute the MCB of the set of possible solutions for uncertain linear systems. Its center can be considered as a candidate. Furthermore, we use the recent results of Zhen et al. (2016) for ARO to characterize for which convex bodies the obtained MCB is optimal.
3. We compare our new method both theoretically and numerically with the RLS method. We show that the solutions x RL S are scale sensitive and may even be outside X . Applications to the input-output model, Colley's Matrix Rankings, and Article Influence Scores demonstrate the advantages and disadvantages of the two methods.
The remainder of the paper is organized as follows. Section 2 discusses the properties of the solution sets, and derive equivalent convex representations of X and X con . The method for computing the MCB is discussed in Sect. 3. In Sect. 4, we compare our method with the RLS method theoretically, and Sect. 5 presents numerical results.
Notation We use [n], n ∈ N to denote the set of running indices, {1, . . . , n}. We use bold faced characters such as x ∈ R n to represent vectors. We use x i to denote the i-th element of the vector x. We denote a · j as the j-th column of the matrix A. We use normal and calligraphy capital letters such as A ∈ R m×n and X to represent matrices and sets, respectively. We denote ζ ∈ R n ζ as the uncertain parameters.

Convexity of united solution set
Let us consider system (5) of uncertain linear equations. We assume that the uncertainty set U that ζ = [ζ 0 ζ 1 · · · ζ n ] ∈ R n ζ resides in, is bounded and defined as follows: and the function f jk is convex in ζ j , j ∈ [n] ∪ {0}, k ∈ [n]. The components of ζ j ∈ U j may be dependent. For i = j, the components of ζ i ∈ U i and ζ j ∈ U j are independent. The uncertainties in the system (5) are indeed column-wise. Note that the dimensions of the uncertain parameters ζ i and ζ j are not necessarily the same for i = j. For a given pair I, J ⊆ [n], where I ∪ J = [n], and I ∩ J = ∅, the corresponding orthant is defined as Lemma 1 The intersection of R n I,J and the solution set X with the uncertainty set U defined in (5): Proof Since this proof is almost identical to the proof for Blanc and Hertog (2008, Proposition 1), we omit it here. Although Blanc and Hertog (2008, Proposition 1) only considers x in the nonnegative orthant with polyhedral uncertainties on the left hand side, their proof also holds for general convex uncertainties on both left and right hand side with x reside in any given orthant. Soyster (1973) shows that the nonnegative sum of a finite number of convex sets is convex, which is basically the result of Lemma 1. Examples 1 and 2 show that in order to preserve convexity of the solution set X the two conditions are indeed necessary: (a) the feasible solutions x ∈ X are within a particular orthant of R n , and (b) the uncertainties in A(ζ ) and b(ζ ) are column-wise.
Example 1 The union of the solutions x ∈ X in two orthants can be nonconvex. Let us consider the following solution set of an uncertain linear system The solution set X can be represented as: One can easily see that the intersection of the solution set X and each orthant of R 2 is indeed convex. However, the set X is nonconvex, see Fig. 1. This result is known in Oettli (1965).
Example 2 The solution set X can be nonconvex when the uncertainties are not column-wise. Let us consider the following solution set of an uncertain linear system: Note that the uncertainties of the system are not column-wise. The set can be represented as: This set is nonempty only in the nonnegative orthant R 2 + , and clearly, the set X is nonconvex. See Fig. 2. A similar observation is also made in, e.g., Tichatschke et al. (1989) and Popova and Krämer (2007).

Convex representation of united solution set
Given the uncertainty sets defined in (8), the intersection of the solution set X and an arbitrary orthant R n I,J can be compactly represented as follows: where the components of the vector functions a · j , b are affine in ζ j , and f jk is convex in ζ j , for all j, k. Due to the presence of products of variables (e.g., ζ j x j for some j ∈ [n]), the representation of set (9) is nonconvex.

Theorem 1
The set (9) admits the following convex representation where 0a · j x j and 0 f jk x j are the recession functions of a · j and f jk , respectively.
Proof The set (10) is obtained by substituting y j = x j ζ j and multiply the inequality constraints containing ζ j with x j in (9). Since the vector functions a · j (·), j ∈ [n], and b(·) are affine, the equalities n j=1 a · j y j x j x j = b(ζ 0 ) are affine in x, ζ 0 and y j , j ∈ [n]. Dacorogna and Maréchal (2008) show that, for a convex function f : R n → R, its perspective g( y, Here is a short proof of the convexity of x f ( y x ) on R n × R + \ {0} uses convex analysis (adopted from Gorissen et al. 2014): from which it follows that g( y, x) is jointly convex since it is the pointwise supremum of functions which are linear in x and y. Hence, the functions x i f ik Theorem 1 shows that for general convex functions f jk , for all j, k, the solution set X in any orthant of R n is convex, which coincides with the findings in Lemma 1.
Our result generalizes the results of Rohn (1981), where the author provides a convex representation of interval linear systems with prescribed column sums.
This transformation technique used in the proof of Theorem 1 is first proposed by Dantzig (1963) to solve Generalized LPs. In the field of disjunctive programming, Balas (1998) employs this technique to derive a convex representation of the union of polytopes. Grossmann and Lee (2003) use it to derive a convex representation of the convex hull of general convex sets. This technique is also applied to the dual of LPs with polyhedral uncertainty in Römer (2010). Gorissen et al. (2014) use it to derive tractable robust counterparts of a linear conic optimization problem. We illustrate this transformation by the following interval linear system example. This example is used throughout this paper.
Example 3 The convex representation of X ∩ R n + with product of variables. Let us consider the solution set in non-negative orthant R 2 [2],∅ = R 2 + : x 1 and ζ 2 = y 2 x 2 , and multiplying the inequality constraints containing ζ j with x j , yields the following representation: One can further simplify this set by eliminating the equality constraints. From Fig. 3, we observe that the set defined in (11) is a full-dimensional polytope.
For interval linear systems,  show that checking the boundedness of the solution set is NP-hard. If we only focus on the solution set in a specific orthant, the boundedness can simply be checked by maximizing i∈I x i − j∈J x j over X in R n I,J .

Convex representations of controllable and tolerable solution sets
Let us first consider a controllable solution set in an arbitrary orthant with column-wise uncertainties:  (12) can be reformulated as: Let us now focus on a special class of (13), where f ik , i, k ∈ [n], are affine. The set (13) can be interpreted as a solution set of a two-stage robust linear optimization problem. Here, the auxiliary variables y j , j ∈ [n], can be seen as general functions of ζ 0 . Ben-Tal et al. (2004) show that solving ARO problems is generally NP-hard, because the auxiliary variables (a.k.a., adjustable variables) are decision rules instead of finite vectors of decision variables. Zhen et al. (2016) show via Fourier-Motzkin elimination that the solution set of two-stage robust linear optimization problems can be reformulated into a convex set, e.g., the set (13) is a polyhedron if U 0 is polyhedral. For computational tractability, we often restrict the adjustable variable y (for simplicity, here we drop the subscript of y j ) to linear decision rules, e.g., the decision rule y is linear in ζ 0 : where the coefficients u ∈ R n y and V ∈ R n y ×n will be optimization variables. Despite that linear functions may not be optimal, it appears that such a decision rule performs well in practice. For instance, Ben-Ameur et al. (2016) and Zhen et al. (2016) show that linear decision rules are optimal for adjustable robust linear optimization problems with a simplicial uncertainty set. Suppose U 0 = {ζ 0 ∈ R n + | 1 ζ 0 ≤ 1} is a standard simplex, then linear decision rules are optimal, and (13) can be reformulated into a polyhedron. For a general convex U 0 , linear decision rules result in an inner approximation of (13).
In case of tolerable solution set, one can eliminate all the auxiliary variables b, and obtain: Note that here we do not restrict the solution set to be within a specific orthant. The set (15) can be seen as a solution set of a static robust optimization problem. If A and B are ellipsoidal, then the set (15) (2012), where the author shows that the tolerable solution set with interval uncertainties is a convex polyhedron.

Maximum size inscribed convex body of solution set
Firstly, in Sect. 3.1, we extend the method of Zhen and den Hertog (2017) to compute the maximum size inscribed convex body (MCB) of a polytopic projection. Since the obtained MCB is an under approximation of the optimal MCB, in Sect. 3.2, we briefly discuss a simple procedure that provides an upper approximation.

MCB of polytopic projection
It is well-known that for a general convex set, finding the MCB can be computationally intractable. In this subsection, we focus on polyhedral solution sets. For instance, the united solution set X ∩ R n I,J defined in (10) is polyhedral if the functions f jk and f 0k are affine, j, k ∈ [n]: for some scalar c jk , c 0k , and vectors d jk , d 0k . Similarly, the controllable and tolerable solution sets can also be reformulated (exactly or approximately via linear decision rules and robust optimization techniques) into polyhedral sets, and the method discussed in this section can also be directly applied. For clarity, we represent the (projected) polyhedral set (16) as follows: where D ∈ R m×(n+n y ) , and c ∈ R m . The auxiliary variable y in (17) represents the variables y j 's and ζ 0 in (16). If the set H is not full-dimensional, then there are (hidden) equality constraints in H. One can simply eliminate the equality constraints via Gaussian elimination. Without loss of generality, let us assume that the set H is full-dimensional.
The set H contains both variables x and y, but it is desired to find the MCB only with respect to x. Due to the existence of y in the description of H, finding the MCB in a polytopic projection is generally a non-convex optimization problem (see Zhen and den Hertog 2017). One can use elimination methods, e.g., Fourier-Motzkin elimination (see Fourier 1824; Motzkin 1936), to eliminate all y in H. This is equivalent to deriving a description of H that does not contain y. Tiwary (2008) shows that deriving an explicit description of a projected polytope is NP-hard. Alternatively, Zhen and den Hertog (2017) approximate the maximum volume inscribe ellipsoid (MVE) of H via linear and quadratic decision rules.
We extend the method of Zhen and den Hertog (2017) to compute the MCB (and its center) of a polytopic projection by solving the following adjustable robust optimization problem: where E ∈ S n with implicit constraint E 0, S n is the set of n × n symmetric matrices, x are non-adjustable variables, the vector function y are called decision rules, and ⊆ R n is a compact convex body. By maximizing the concave function logdet(E) in Problem (18), the matrix E "stretches" (or "suppresses") and "rotates" the convex body around x in H to its maximum volume. If one would like to "stretch" (or "suppress") the convex body along the x's axes (i.e., without rotation), this can simply be done by restricting the non-diagonal entries of E to be 0, and maximizing i∈ [n] log E ii . When H is unbounded, the volume of the MCB is also unbounded. The boundedness of H can be checked via an LP problem.
Let us restrict the decision rule y to be linear in ξ as in (14), and substitute it in (18) to get the affinely adjustable robust formulation: Problem (19) where ∈ R m×m ξ . Let us consider two special polyhedral sets, i.e., is a box or simplex. Zhen et al. (2016) proves that linear and two-piecewise decision rules are optimal for two-stage robust linear optimization problems with simplex and box uncertainties, respectively. If = {ξ ∈ R n + | 1 ξ ≤ 1} is a standard simplex, the solution of (20) gives a maximum inscribed simplex of H; if = {ξ ∈ R n | |ξ | ≤ 1} is a box, there exists two-piecewise affine decision rules that are optimal for (18), and the solutions of (20) are in general suboptimal. For tolerable solution sets with interval uncertainties, Hladík and Popova (2015) devise an exponential LP methods for computing maximal inner boxes exactly, and also propose a polynomial heuristic that yields an approximated maximal inner box. Furthermore, Ben-Tal et al. (2004) show that static decision rules are optimal for two-stage robust linear optimization problems with constraint-wise uncertainties. Hence, one can directly show that the solution of where d i· ∈ R n is the i-th row of the matrix D. Problem (21) approximates the MVE of H. A largest ball in H can be determined by replacing the matrix E by eI ∈ R n×n (i.e., a product of a scaler e ∈ R + and identity matrix) everywhere in (21). For the rest of this paper, we focus on MVE, i.e., = {ξ ∈ R n | ||ξ || 2 ≤ 1}, because it possesses many appealing properties: (a) it is unique, invariant of the representation of the given convex body; (b) its center is a centralized (relative) interior point of the convex body; (c) it is a inner approximation of H which admits a simple description (i.e., one inequality constraint). We denote x aMV E as the approximated MVE center of H obtained from the MVE method (21). In the following example, we solve (21) to compute x aMV E of the solution set H.
Example 4 The maximum volume inscribed ellipsoid (Example 3 continued). From Fig. 3 we know that H is full-dimensional. Since we are interested in the MVE center of H only with respect to x, we compute the x aMV E of H by solving (21), and obtain x aMV E = (52.1, 30.7). In order to evaluate the obtained solution, we derive an explicit description of H with no auxiliary variables by using Fourier-Motzkin elimination, and obtain the optimal MVE center x MV E at (53.6,30). Note that, in general, it is NP-hard to derive such a description. From Fig. 3, one can observe that x aMV E is a close approximation of x MV E . Hadjiyiannis et al. (2011) propose to compute upper bounds on the optimal value of adjustable robust optimization problems by only considering a finite set of scenarios, which they call critical scenarios (CSs). The CSs are obtained by solving the auxiliary optimization problems:

Upper bounding method of Hadjiyiannis et al. (2011)
where E * and V * denote the optimal solution from (19). If more than one CS is determined from the k-th constraint, an arbitrary CS is chosen and included in the CS set. The scenario counterpart of Problem (18) with respect to the CS set , where = {ξ 1 , . . . , ξ m }, is given by the following optimization problem: For each ξ k ∈ , k ∈ [m], we only need a feasible y k ∈ R n y to exist. Since ⊂ , Problem (23) provides an upper bound of (18).

Comparison of solution methods
In this section, we compare the theoretical aspects of the robust least-square (RLS) method with our MVE method discussed in Sect. 3. Let us first consider an interval linear system where: and ζ = [ζ 0 ζ 1 · · · ζ n ] ∈ R n ζ . We further denote the nominal realization ζ nom as the median of the intervals, i.e., ζ nom = 1 2 (ζ + ζ ) ∈ R n ζ . (24), the solution x RL S is optimal for (7) if and only if it is also optimal for the following SOCP problem with a unique z ∈ R m :
Note that in Lemma 2, the optimal solution x RL S of (25) is not restricted to be within any specific orthant. The tractability of Problem (7) strongly relies on the choice of the uncertainty set U. For deriving tractable counterpart of (7) with ellipsoidal uncertainties in A and/or b, we would like to refer to Ghaoui and Lebret (1997), Beck and Eldar (2006) and Jeyakumar and Li (2014). In the following example, we solve Problem (25) for the interval linear system in Example 3.
Example 5 The robust least-squares solution for an interval linear system. We apply the RLS method to the interval linear system in Example 3 and find the solution x RL S is at (67.06, 10.59), which is denoted as in Fig. 3. The solution x RL S coincides with the nominal solution x nom (i.e., the solution obtained from the nominal realization) of the system.
The RLS method is in line with the philosophy of Robust Optimization (see Ben-Tal et al. 2009), i.e., minimizing the violation with respect to the worst-case scenario. Our method finds a centered solution of the solution set. In Sect. 5, we show that in many real-life problems that involve solving a system of linear equations, the uncertainties are often column-wise, i.e., the uncertainties within each column of A(ζ ) and b(ζ ) may be dependent, and the centered solution x aMV E can be efficiently obtained. However, the choice of uncertainty sets for Problem (7) that admits a tractable counterpart is rather limited, e.g., for polyhedral uncertainties, Problem (7) is generally NP-hard.
One of the most fundamental properties of a system of (uncertain) linear equations is scale invariance. The nominal solution x nom and centered solution x aMV E are scale invariant, whereas, the solution x RL S is not. In fact, x RL S can be outside the solution set (see Example 6), and even if one restrict x RL S to be within X , it may still be on the boundary of the solution set. As the feasibility of the solutions are not guaranteed, the solution x RL S may be less appealing than x aMV E .
Example 6 Scale sensitivity of the solutions. Let us consider an adapted version of the interval linear system in Example 3 where the uncertainty sets are defined as: The components of the first row of the interval linear system in Example 3 are now multiplied by a factor 30. Note that this operation does not alter the set X ∩ R 2 + . The solutions for this uncertain linear system are depicted in Fig. 4. The nominal and MVE solutions remain unchanged. The solution x RL S is outside the solution set.

Numerical experiments
In this section, we conduct four experiments to evaluate the candidate solutions. The first experiment considers the interval linear system introduced in Example 3 and its adapted version in Example 6. The other three are input-output model, Colley's Matrix Ranking and Article Influence Scores, respectively. A common feature of these problems is that their solutions are obtained by solving a system of linear equations. All the procedures are performed by using Mosek 8 (see MOSEK ApS 2017) within Matlab R2017a on an Intel Core i5-4590 CPU running at 3.3 GHz with 8 GB RAM under Windows 7 operating system.

A simple experiment
Firstly, we compute the candidate solutions for the interval linear system discussed in Example 3. Given a candidate solutionx, three measures are considered: , are obtained from the Hit-and-Run sampling (see Smith 1984).
For the interval linear system in Example 3, the nominal solution x nom coincides with x RL S (see Fig. 3). From Table 1, one can observe that x nom and x RL S are most robust against the worst-case deviations. The solutions x MV E 1 and x aMV E are centered solutions (see also the exact ranges of x in the last column of Table 1). The In Table 2, we evaluate the solutions for the interval linear system in Example 6. Here, the nominal solution x nom is no longer the same as x RL S , and x RL S is outside the solution set (see Fig. 4). Therefore, the corresponding volume of the MVE is 0. The solutions x nom , x MV E and x aMV E are scale invariant.

Production vector for input-output model
Leontief's Nobel prize-winning input-output model describes a simplified view of an economy. Its goal is to predict the proper level of production for each of several types of goods or service. We apply this to predict the production of different industries in the Netherlands. In Table 3, we present the data that are reported in Deloitte (2014). This is a simplified version of the consumption data of the Netherlands published by the Dutch statistics office. From Leontief (1986), the nominal input-output matrix is defined as: where w ∈ R 5 is the total output vector, C ∈ R 5×5 is the consumption matrix from Table 3, and Diag(·) places its vector components into a diagonal matrix. The nominal production vector x nom can be obtained by solving the following system of linear equations: where b is the vector of the nominal external demands (see the last column of Table 3). Suppose there are uncertainties in the system (26), and each component of A(ζ ) = [ζ 1 · · · ζ n ] and b(ζ ) = ζ 0 resides in an independent interval, where ζ = [ζ 0 ζ 1 · · · ζ n ] ∈ R n 2 +n . We assume that the interval uncertainty sets U j , j = [n] ∪ {0}, are as follows: Five industries are considered, i.e., agriculture, fishing, forestry (AFF) industry, manufacturing industry, service industry, education and healthcare (E & H) industry and other industries. The external demand (ED) and the total output are also reported Table 4 The candidate production vectors for σ = 15% x nom The solutions from the MVE method and the RLS method are denoted as x aMV E and x RL S , respectively.
The solution x nom is the nominal solution. The exact ranges of the components of x are reported in the last column. LS denotes the complexity of solving a linear system where σ is user specified, a · j is the j-th column of the nominal matrix A. Inputoutput models with interval data were also studied in, e.g., Rohn (1978) and Dymova et al. (2013). In Table 4, the computation time of the solution methods is positively correlated with its theoretical complexity. Since the problem size is relative small, all the solutions can be obtained within 2 s. We again consider the three measures as in Sect. 5.1. From Table 5, we observe that the solution x aMV E has the largest inscribed ellipsoid and the best (i.e., least) M D; x RL S is the most robust solution with respect to W D, but geometrically, it is not centered. The nominal solution x nom is the second best solution for all the three measures. For other values of σ , or different σ i j 's for each components of A and b, the above observations remain valid. However, when σ ≥ 25%, the RLS solutions are outside the solution set. Since x aMV E is an approximation of the optimal MVE center, we apply the upper bounding method of Sect. 3.2. The upper bounding V O L is 44.4. The obtained inner approximation is 44.3, which implies that x aMV E is very close to the optimal solution. The solutions from the MVE method and the RLS method are denoted as x aMV E and x RL S , respectively.
The solution x nom is the nominal solution. The bold numbers show that the corresponding solution performs the best (among the candidate solutions) with respect to the corresponding measure

Colley's matrix ranking
Colley's bias free college football ranking method was first introduced by Colley (2001). This method became so successful that it is now one of the six computer rankings incorporated in the Bowl Championship Series method of ranking National Collegiate Athletic Association college football teams. The notation here is adapted from Burer (2012).
Colley Matrix Rankings require to solve a system of linear equations Ax = b. For n teams, the n × n matrix W is defined as W i j = number of times team i has beaten team j.
In particular, W i j = W ji = 0 if i has not played against j, and W ii = 0 for all i. Note that the i j-th entry of W + W represents the number of times team i and team j has played against each other. Let 1 be the all-ones vector, then the i-th entry of (W + W )1 and (W − W )1 gives the total number of games played by team i, i.e., the schedule of the games, and its win-loss spread. The Colley matrix A and the vector b are defined via the schedule of the games and the win-loss spread vector respectively, i.e., where I is the identity matrix and Diag(·) places its vector components into a diagonal matrix. Since the schedule of the games is often predetermined, we only consider uncertainties in the vector b. We empirically investigate the robust version of Colley Matrix ratings to modest changes in the win-loss outcomes of inconsequential games. A game is inconsequential if it has occurred between two bottom teams, i.e., teams win less than 50% of all the games they played. Suppose m inconsequential games has been played during the whole season. Let ζ ∈ R m denote the perturbation of the games. The game j switches its outcome if ζ j = 1, and it remains unchanged if ζ j = 0. For all j, we have 0 ≤ ζ j ≤ 1. Then, we define a matrix ∈ R n×m , where i j = ⎧ ⎨ ⎩ 1 if team i loses the game j −1 if team i wins the game j 0 otherwise. The vector ζ represents the possible switches in the outcome of the games. The maximum number of inconsequential games that are allowed to switch their outcomes is less than L ∈ N 0 , i.e., j ζ j ≤ L. The polyhedral solution set is as follows: where the matrix A and the vector b are nominal, conv(X ) denotes the convex hull of the set X , and U = ζ : ζ ∈ {0, 1} m , j ζ j ≤ L . The uncertainty set U contains all possible integral ζ 's (i.e., scenarios). For ζ = 0, the nominal rating vector x nom = A −1 b is on the boundary of X . Note that the ratings of Colley's Matrix Rankings are not necessarily nonnegative. Since negative ratings are rather rare and their values are marginal (often very close to zero), we restrict ourselves to the rating vectors that are nonnegative.
The data we use in this subsection are downloaded from the website Wolfe (2017). The data contains the outcomes of 4197 college football games played by n = 760 teams in 2016, US. There are m = 112 inconsequential games in total. We allow at most L = 30 inconsequential games to switch their outcomes. The total number of possible outcomes is |U| = 30 i=1 112 i = 2.4 × 10 27 . The solution x aMV E is defined as the approximated MVE center of the solution set conv(X ). For the polyhedral set conv(X ), the RLS method requires solving a 2-norm maximization problem which is NP-hard. Burer (2012) proposes a two-stage method to solve the following MINLP problem: We denote the solution of Burer as x RL S . Due to the high dimension of the solutions (i.e., n = 760), we do not report the obtained solutions, and the exact ranges of the components of x for this numerical experiment. Since the uncertainty set U is discrete, the measures that we consider here are slightly different from those introduced in Sect. 5.1: -V O L the volume of the MVE centered atx within conv(X ) -W D the worst-case 2-norm deviations of Ax from b + ζ with respect to 10 5 randomly sampled integer ζ ∈ U (i.e., max ζ ∈U ||Ax − b − ζ || 2 ) -M D the mean 2-norm deviations of 10 5 uniformly sampled solutions in conv(X ) (i.e., 1 n s i∈[n s ] ||x i −x|| 2 ), where x i ∈ X , i ∈ [n s ], are obtained from the Hit-and-Run sampling (see Smith 1984) -M D D the mean 2-norm deviations of Ax from b + ζ with respect to 10 5 randomly sampled discrete ζ ∈ U.
From Table 6, it is readily obvious that the solution x aMV E is the best for three out of four measures. Due to the problem definition, the effect of switching the result of the inconsequential games is not symmetric. The solution x nom is on the boundary of X and it is the least robust solution among all three solutions with respect to the considered measures. We again evaluate the quality of the approximation x aMV E by computing its upper bounding volume. The obtained upper bounding volume is 0.18. The optimal volume lies between 0.06 and 0.18. The MINLP problem is solved with CPLEX 12.7 ILOG (2016). For larger sized problems, one can expect exponential growth in computation time for the RLS method, whereas the MVE method remains computationally tractable.

Article influence scores
Around 1996-1998, Larry Page and Sergey Brin, Ph.D. students at Stanford University, developed the PageRank algorithm for rating and ranking the importance of Web pages (see Brin and Page 1999). An adapted version of PageRank has recently been proposed to rank the importance of scientific journals as a replacement for the traditional impact factor (see Bergstrom et al. 2008 (28) Finally, we construct the matrix A, a convex combination of S and a rank-one matrix, i.e., where α is the damping factor and w1 is a n × n matrix. The damping factor models the possibility that a searcher choose a random paper out of all papers. Therefore, the closer the α gets to 1, the better the journal's citation structure is represented by the matrix A. The influence vector x * can be obtained by solving as follows: From the Perron-Frobenius theorem, we know a unique rating vector x * can be found. The Article Influence score of journal i can be calculated as follows: ∀i.
In this subsection, we assume α = 90%. Let us consider the matrix S and vector w defined in (28) and denote the obtained matrix in (29) as the nominal matrix A. The nominal influence vector x nom is obtained by solving the system of linear equations (30). Since the estimated probabilities are not exact, we take uncertainty in the matrix A into consideration. Let us consider the following column-wise 1-norm uncertainties in A: U = ζ : ||ζ j − a · j || 1 ≤ σ, ζ j 1 = 1, ζ j ≥ 0, ∀ j , where ζ = [ζ 1 · · · ζ n ] ∈ R n 2 , σ = 20%, ζ j and a · j are the jth column of matrix A(ζ ) and A, respectively. The uncertainties occur in the left-hand side of the system. Note that each column of the nonnegative matrix A(ζ ) is a probability vector, i.e., ζ j 1 = 1 for all j. Hence, the uncertain parameters are dependent. Since 2-norm The x nom denotes the nominal solution; the x aMV E is the approximated MVE center obtained by solving (21). The bold numbers show that the corresponding solution performs the best (among the candidate solutions) with respect to the corresponding measure. LS denotes the complexity of solving a linear system maximization over a polyhedron is an N P-hard problem, the RLS method is computationally intractable. We consider the tractable (upper-bound) approximation of the RLS proposed in Juditsky and Polyak (2012): The solution of this approximation coincides with x nom . We refer to Polyak and Timonina (2011) for a fast algorithm that solves (J P) with high-dimensional A. Hence, in the remaining of this section, we do not distinguish the solution of (J P) from x nom . The resulting Article Influence Scores from the influence vectors are reported in Table 7. The exact ranges of the components of the solution and the AIS AI are reported. The difference between the nominal and MVE solution is marginal. The width of [x, x] and [AI , AI ] indicates that the system (30) is sensitive to this type of uncertainties. We again consider the measures V O L and M D. Besides these two measures, the mean 2-norm deviations of 10 4 uniformly sampled (A, b) in U are also considered (i.e., M D A,b ). From Table 8, one can observe that the solution x aMV E is slightly more robust than x nom with respect to all three considered measures. Here, the nominal solution is robust against uncertainties. The obtained upper bounding volume of x aMV E is 0.0739. The optimal MVE volume lies between 0.0313 and 0.0739. We further observe that for a smaller uncertainty σ or damping factor α, the difference between x nom and x aMV E is smaller.

Conclusion and future research
In this paper, we first generalize the results for interval linear systems. For a system of uncertain linear equations with column-wise uncertainties, we derive convex representations of the united solution set in a given orthant. The exact ranges of the components of the solutions can then be determined. Via a convex representation technique and the techniques from (adjustable) robust optimization, we show how to derive convex representations of a broad class of controllable and tolerable solution sets. We apply the MCB method for obtaining centered solutions of systems of uncertain linear equations, and compare our proposed method both theoretically and numerically with the RLS method. The solutions from the RLS method may even be outside the solution set. As a byproduct, the MCB method produce a simple inner approximation of the solution set. From the numerical experiments, we observe that, for column-wise dependent uncertainties, our proposed solutions are more centered than the RLS or nominal solutions.
Further research is needed to determine the usefulness the MCB method for many other real-life applications, e.g., analysis of mechanical structures, electrical circuit designs and chemical engineering.