- Email: [email protected]

Contents lists available at SciVerse ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Numerical range of some doubly stochastic matrices q Kristin A. Camenga a, Patrick X. Rault b, Daniel J. Rossi b, Tsvetanka Sendova c, Ilya M. Spitkovsky d,⇑ a

Department of Mathematics and Computer Science, Houghton College, Houghton, NY 14744, USA Department of Mathematics, State University of New York at Geneseo, Geneseo, NY 14454, USA c Department of Mathematics and Computer Science, Bennett College, Greensboro, NC 27401, USA d Department of Mathematics, College of William and Mary, Williamsburg, VA 23187, USA b

a r t i c l e

i n f o

a b s t r a c t A classiﬁcation of all possible shapes is given for numerical ranges of 4 4 doubly stochastic matrices. The tests determining the shape are also provided, along with illustrating examples. Ó 2013 Elsevier Inc. All rights reserved.

Keywords: Numerical range Doubly stochastic matrices

1. Introduction With F standing for either the ﬁeld R of real or C of complex numbers, we denote by Fnm the linear space (algebra, if m ¼ n) of the n m matrices with entries in F. We also use the standard abbreviation Fn :¼ Fn1 , and supply Cn with the inner P product hx; yi ¼ nj¼1 xj yj and the respective norm kxk :¼ hx; xi1=2 . The numerical range of A 2 Cnn by deﬁnition is

WðAÞ ¼ fhAx; xi : kxk ¼ 1g: This is a compact convex subset of C containing the spectrum rðAÞ of A; see e.g. [1] or [2, Chapter 1] for these and other fundamental properties of WðAÞ. We will repeatedly be using the following: (i) WðAÞ is unitarily invariant: WðU T AUÞ ¼ WðAÞ for any unitary U 2 Cnn . (ii) If A is unitarily reducible, that is, unitarily similar to the direct sum A1 A2 , then WðAÞ is the convex hull of WðA1 Þ [ WðA2 Þ; that is WðAÞ ¼ convfWðA1 Þ; WðA2 Þg. From these properties it follows in particular that WðAÞ ¼ conv rðAÞ whenever A is normal. On the other hand, WðAÞ is an ellipse with the foci at the eigenvalues of A if A 2 C22 and is not normal (the Elliptical Range Theorem). The case n ¼ 3 has also been treated fully, see [3] and its English translation [4] for the classiﬁcation, and [5–7] for the pertinent tests. However, starting with n ¼ 4, the shape of WðAÞ is known only for some special classes of matrices. Here we concentrate our attention on doubly stochastic matrices, that is, matrices A 2 Rnn for which the entries are nonnegative while the row and column sums are all equal to one: n X aij ¼ 1 for i ¼ 1; . . . ; n; j¼1

q

and

n X aij ¼ 1 for j ¼ 1; . . . ; n:

ð1:1Þ

i¼1

The work on this paper was begun during the REUF workshop at the American Institute of Mathematics in July 2011, supported by the NSF.

⇑ Corresponding author.

E-mail addresses: [email protected] (K.A. Camenga), [email protected] (P.X. Rault), [email protected] (T. Sendova), [email protected] wm.edu, [email protected] (I.M. Spitkovsky). 0096-3003/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2013.06.011

41

K.A. Camenga et al. / Applied Mathematics and Computation 221 (2013) 40–47

Observe that (1.1) can be interpreted as e ¼ ½1; . . . ; 1T being an eigenvector of both A and AT , corresponding to the eigenvalue one. Thus, doubly stochastic matrices are unitarily reducible to the direct sum of ð1Þ with some A1 2 Rðn1Þðn1Þ :

U T AU ¼ ð1Þ A1 ;

ð1:2Þ

where U is an orthogonal matrix with the ﬁrst column equal p1ﬃﬃn e. Consequently, WðAÞ in this case is the convex hull of WðA1 Þ with the point 1. Moreover, 1 is a corner point of WðAÞ [8]. If A (equivalently, A1 ) is normal, then WðAÞ is the convex hull of rðAÞ ¼ rðA1 Þ [ f1g. Thus we are only interested in the case of non-normal A1 . This, combined with the Elliptical Range Theorem, immediately implies a complete description of numerical ranges for 3 3 doubly stochastic matrices, normal or otherwise; see [9] for details. Our aim in this paper is to realize a similar plan for n = 4. According to [4], for a non-normal A1 there are a priori three possibilities: (i) WðA1 Þ is the convex hull of a point and an ellipse (with the point lying either inside or outside the ellipse), (ii) the boundary of WðA1 Þ contains a ﬂat portion, with the rest of it lying on a 4th degree algebraic curve, or (iii) WðA1 Þ has an ovular shape, bounded by a 6th degree algebraic curve. Since WðA1 Þ also is symmetric with respect to the real axis (as A1 has real entries), the ellipse in case (i) has one of its axes on R, and the ﬂat portion in case (ii) is vertical. In the forthcoming sections we provide criteria to distinguish between the three possibilities, show that each of them actually materializes, and describe the respective shapes of each WðAÞ. To simplify the terminology, we label the shape of WðAÞ by the degree of a non-ﬂat portion of its boundary. 2. Numerical range of 2nd degree In this section we characterize the case when the non-ﬂat portion of @WðAÞ is an elliptical arc. Theorem 1. Let A be a 4 4 doubly stochastic matrix. Then @WðAÞ consists of elliptical arcs and line segments if and only if

l :¼ tr A 1 þ

tr A3 trðAT A2 Þ

ð2:1Þ

trðAT AÞ tr A2

is an eigenvalue of A (multiple, if

tr A 1 3l > 0;

l ¼ 1). If, in addition,

ðtr A 1 3lÞ2 tr ðAT AÞ þ 1 þ 2ðdet AÞ=l þ l2 > 0;

ð2:2Þ

then WðAÞ also has corner points at l and 1, and thus four ﬂat portions on the boundary. Otherwise, 1 is the only corner point of WðAÞ, and @WðAÞ consists of two ﬂat portions and one elliptical arc. tr A3 trðAT A2 Þ

1 1 Proof. Due to (1.2), the right hand side of (2.1) coincides with tr A1 þ tr ðA1T A Þtr . From [7, Theorem 5], it then follows that A21 1 1 WðA1 Þ is the convex hull of an ellipse and a point if and only if l is an eigenvalue of A1 . This proves the ﬁrst statement of the theorem. Further, recall from [5] that if the above condition on l is satisﬁed, then WðA1 Þ ¼ convfE; lg, where E is an ellipse with foci at two other eigenvalues k1 ; k2 of A1 and the minor axis d ¼ tr AT1 A1 jk1 j2 jk2 j2 l2 . So, WðAÞ ¼ convfE; l; 1g has two corner points (at 1 and l) if l lies to the left of E, and only one corner point (at 1) otherwise. Now observe that l lies to the left of E if and only if it is located in the exterior of E, that is,

ðjk1 lj þ jk2 ljÞ2 jk1 k2 j2 > d;

ð2:3Þ

and to the left of its center ðk1 þ k2 Þ=2. The second condition is equivalent to

3l < k1 þ k2 þ l; where the right hand side is of course tr A 1. Collecting all the terms on the same side of the inequality, we arrive at the ﬁrst inequality in (2.2). It remains to show the equivalence of its second inequality and (2.3). To this end, we consider separately the case of complex conjugate and real kj . Case 1. k1 ¼ k2 ¼: a þ bi. Inequality (2.3) then can be rewritten as

2 2 2 4 ða lÞ2 þ b ð2bÞ > trðAT1 A1 Þ 2ða2 þ b Þ l2 :

Simplifying the left side and right side individually yields

ðtr A1 3lÞ2 > trðAT1 A1 Þ 2ðdet A1 Þ=l l2 : Putting this in terms of A, we arrive at the desired second inequality in (2.2). Case 2. All eigenvalues of A1 are real. Recalling that l lies to the left of the eigenvalues, inequality (2.3) becomes

ðk1 þ k2 2lÞ2 ðk1 k2 Þ2 > trðAT1 A1 Þ k21 k22 l2 :

42

K.A. Camenga et al. / Applied Mathematics and Computation 221 (2013) 40–47

Working with the left hand side gives

4l2 þ 4k1 k2 4ðk1 þ k2 Þl > trðAT1 A1 Þ k21 k22 l2 : Rewriting this and collecting terms yields

5l2 þ 2k1 k2 4ðtr A1 lÞl þ ðk1 þ k2 Þ2 trðAT1 A1 Þ > 0: After putting this in terms of A and factoring gives

5l2 þ 2ðdet AÞ=l þ ðtr A 1 lÞðtr A 1 5lÞ trðAT AÞ þ 1 > 0: Noting that

ððtr A 1 3lÞ 2lÞððtr A 1 3lÞ þ 2lÞ ¼ ðtr A 1 3lÞ2 4l2 completes the proof.

h

The following examples demonstrate that each case considered by Theorem 1 actually occurs. Example 1. Consider the doubly stochastic matrix

0 B B A¼B @

0

1 3

1 3 1 4 5 12

0 9 32 37 96

1 4 1 2 1 4

0

5 1 12 1 C 6 C : 7 C A 32 19 96

Using equation (2.1), we compute

l ¼ 1=3. The characteristic polynomial of A is

1=128 ð11xÞ=192 ð193x2 Þ=384 ð43x3 Þ=96 þ x4 ; which evaluated at 1=3 is 0. In addition, we compute the formulas in inequalities (2.2) as

tr A 1 3l ¼ 43=96 and ðtr A 1 3lÞ2 trðAT AÞ þ 1 þ 2ðdet AÞ=l þ l2 ¼ 779=9216: Since the latter is negative, Theorem 1 states that @WðAÞ has two ﬂat portions and one elliptical arc. Indeed, A reduces unitarily to ð1Þ A1 for some 3 3 matrix A1 , and in Fig. 1 we give WðAÞ, and the horizontal ellipse it contains, WðA1 Þ. Example 2. Next, we give a similar example where the elliptical arc is part of a vertical ellipse:

0 B B A¼B @

0

1 3

1 3 1 4 5 12

0 1 8 13 24

1 4 1 2 1 4

0

5 1 12 1 C 6 C : 3 C A 8 1 24

By coincidence, we again compute

l ¼ 1=3, and though the characteristic polynomial is now

2

ð1=32Þ ð3xÞ=16 ð47x Þ=96 ð7x3 Þ=24 þ x4 ; it again has l as a root. The formulas in inequalities (2.2) evaluate to 7=24 and 59=576 respectively, so the number of ﬂat portions is still the same. However, Fig. 2 shows that WðAÞ is the convex hull of a vertical ellipse and the point 1, as opposed to the horizontal ellipse in the previous example. Indeed, the eigenvalues marked in the graph include a complex conjugate pair.

Fig. 1. Numerical range of A with two ﬂat portions and an elliptical arc on the boundary, and the ellipse it contains. The eigenvalues of A are indicated by the points.

K.A. Camenga et al. / Applied Mathematics and Computation 221 (2013) 40–47

43

Fig. 2. Numerical range A with two ﬂat portions and an elliptical arc on the boundary, and the contained vertical ellipse. The eigenvalues of A are indicated by the points.

Fig. 3. Numerical range of A with four ﬂat portions and elliptical arcs on the boundary, and the numerical range of the unitary reduction A1 , which is the convex hull of a vertical ellipse and l. The eigenvalues of A are indicated by the points.

Example 3. Next, we give an example of a doubly stochastic matrix, A, where @WðAÞ consists of four ﬂat portions and two elliptical arcs. Let a be the second largest root of the four real roots of the polynomial 28995327 þ 108880894x þ 143424512x2 778960896x3 þ 603979776x4 and let b be the second smallest of the four real roots of 5072 þ 17953x þ 310274x2 1013248x3 þ 786432x4 . With

0

0

a

B 1 B A¼B 3 @ b

0

2 3

b

1 256 255 256

1 4 1 2 1 4

a 0

3 4

1

a

1 6 191 b 256 509 768 þbþ

C C C; A

a

the characteristic polynomial evaluated at l is identically 0. The two formulas in inequality (2.2) come out to two positive numbers, roughly 0:65 and 0:31 respectively. Hence WðAÞ is the convex hull of an ellipse and two points outside of the ellipse. Fig. 3 shows WðAÞ, and that the ellipse is vertical. Example 4. Lastly, we give an example of a doubly stochastic matrix, A, where WðAÞ is the convex hull of a horizontal ellipse and two points. Let a be the second smallest root of the polynomial 301228513 þ 853379010xþ 4421068800x2 17180393472x3 þ 14495514624x4 and let b be the second largest root of the polynomial 10016 19031x þ 1045714x2 3505152x3 þ 3145728x4 . Just as in the previous examples, the characteristic polynomial of the matrix

0

0

B 1 B A¼B 3 @ b 2 b 3

a 0 109 512 403 512

1 4 1 2 1 4

a 0

3 4

1

a

1 6 275 b 512 697 1536 þ b þ

C C C A

a

44

K.A. Camenga et al. / Applied Mathematics and Computation 221 (2013) 40–47

Fig. 4. Numerical range of A with four ﬂat portions and elliptical arcs on the boundary, and the numerical range of the unitary reduction A1 , which is the convex hull of a horizontal ellipse and l. The eigenvalues are indicated by the points.

vanishes at l. In this case the two formulas in inequality (2.2) come out to two positive numbers of about 0:77 and 0:21 respectively. Hence WðAÞ is the convex hull of an ellipse and two points outside of the ellipse and has four ﬂat portions. Fig. 4 shows WðAÞ, and that in this case the ellipse is horizontal.

3. Numerical range of 4th degree By Birkhoff’s Theorem, the doubly stochastic matrices form a convex polytope in Rnn , its vertices being the n! permutation matrices. So, it is not surprising that the set of all doubly stochastic matrices A with the ﬁxed symmetric part H ¼ ðA þ AT Þ=2 also forms a polytope. The latter is stated more rigorously in the following lemma. Lemma 2. Consider A 2 Rnn with given symmetric part H:

A ¼ H þ K;

where K ¼ ðA AT Þ=2:

ð3:1Þ

Then A is doubly stochastic if and only if H is doubly stochastic, the row (and thus, column) sums of K are zero, and its upper diagonal entries satisfy

jkij j 6 hij ;

1 6 i < j 6 n:

ð3:2Þ

Proof. Directly from (1.1) it follows that if A is doubly stochastic, then its symmetric part H is doubly stochastic as well. On the other hand, if H is doubly stochastic, then for A given by (3.1) to be doubly stochastic it is necessary and sufﬁcient that K has zero row/column sums and hij kij P 0. The latter condition is equivalent to (3.2). h Note that antisymmetric matrices with zero row/column sums form a ðn1Þðn2Þ -dimensional subspace of Rnn . 2 So, (3.2) can be rewritten as a system of nðn 1Þ=2 linear inequalities in ðn 1Þðn 2Þ=2 variables. In particular, for n ¼ 4 we have

2

0

0

1

1

6 0 0 6 K ¼ x6 4 1 0

0 0

0 1

1

0 1

0

3

2

0

1

7 6 1 0 7 6 7 þ y6 5 4 0 1 1

0

0

1

1 0

0 1

1

0

3

2

0

7 6 1 7 6 7 þ z6 5 4 0 1

1

0 1

3

0 0

0 0

1 0

7 7 7; 5

1 0

0

ð3:3Þ

and (3.2) takes the form

jxj 6 h13 ; jyj 6 h23 ; jzj 6 h24 ; jx þ yj 6 h34 ; jy þ zj 6 h12 ; jx þ y þ zj 6 h14 :

ð3:4Þ

For the following theorem, recall our notation e for the vector with all the entries equal to one. Theorem 3. Let A 2 R44 be a doubly stochastic matrix. Then @WðAÞ consists of 4th degree algebraic arcs and line segments if and only if in its representation (3.1),

H ¼ kI þ

1k T T ee ff 4

ð3:5Þ

for some f 2 R4 orthogonal to e, and k satisfying

1 ð1 þ 4 minffi2 gÞ 6 k 6 1 þ 4minffi fj g; i–j 3

ð3:6Þ

K.A. Camenga et al. / Applied Mathematics and Computation 221 (2013) 40–47

45

?

while K is given by (3.3), satisﬁes (3.4), and has no eigenvectors lying in ff ; eg . The upper signs in (3.5) and (3.6) correspond to the case when @WðAÞ has three ﬂat portions, one of them being vertical and lying to the left of WðAÞ. Alternatively, if (3.5) and (3.6) hold with the lower signs, @WðAÞ has only two ﬂat portions, with 1 being their common endpoint. Proof. In order for @WðAÞ to contain algebraic arc(s) of 4th degree, the same should be true for @WðA1 Þ, with A1 as in (1.2). According to [5, Section 3], this happens if and only if A1 is unitarily irreducible and its real part H1 has a multiple maximal or minimal eigenvalue k. Equivalently, k is a multiple eigenvalue of H from (3.1), which is either minimal or second in value after one. Along with the double stochasticity of H (see Lemma 2), this implies (3.5) and (3.6). Conditions (3.3) and (3.4) on K follow from the description of all doubly stochastic matrices with ﬁxed symmetric part (Lemma 2 again). Finally, A1 is unitarily irreducible if and only if e is the only common eigenvector of H and K, that is, if K has no eigenvectors orthogonal to e and f . Suppose now that all the conditions (3.3)–(3.6) are satisﬁed. Then A1 is unitarily irreducible and, according to [5], its numerical range has a ﬂat portion on its boundary lying on the vertical line x ¼ k. Just as A reduces to ð1Þ A1 , we reduce H to ð1Þ H1 . If k is the minimal eigenvalue of H1 (that is, (3.5) and (3.6) hold with the upper signs), then this vertical portion is to the left of WðA1 Þ and thus of WðAÞ. Therefore, it is part of @WðAÞ, along with two other ﬂat portions, originating at 1. On the other hand, for (3.5) and (3.6) holding with the lower signs, the vertical portion of @WðA1 Þ is to the right of WðA1 Þ and is thus absorbed by the interior of WðAÞ. h To demonstrate that such matrices actually exist, we give examples of numerical ranges with ﬂat portions on the boundary to illustrate the above theorem. Example 5. Consider the following matrix.

02

2 5 1 5

B 0 B A¼B @ 0 35 0

1 5 2 5 2 5

1 5

2 5

0

5 2 5

0

2 5

1 C C C: A

Eq. (3.1) becomes A ¼ H þ K, where H ¼ ðA þ AT Þ=2 and K ¼ ðA AT Þ=2. H has three eigenspaces: the ﬁrst generated by e, with eigenvalue 1, the second generated by ½3; 1; 1; 1T with eigenvalue k2 ¼ 15, and the third a 2-dimensional space gener3 ated by ½0; 1; 0; 1T and ½0; 1; 1; 0T corresponding to the double eigenvalue k3 ¼ 2 ¼ k4 . Computing H k3 I 1k eeT yields 4 5 a symmetric positive semi deﬁnite matrix

0 M¼

9 20 B 3 B 20 B 3 @ 20 3 20

3 3 3 20 20 20 1 20 1 20 1 20

1 20 1 20 1 20

1 20 1 20 1 20

1 C C C A T

of rank 1, and the vector f from its representation M ¼ ff is determined up to a sign and can be found via dividing any of the columns of M by the square root of the respective diagonal entry. Consequently, f ¼ 2p1 ﬃﬃ5 ½3; 1; 1; 1T , and T 9 3 H ¼ k3 I þ 1k eeT þ ff . Furthermore, note that inequality (3.6) is satisﬁed: 81 6 2 6 10 . We conclude that @WðAÞ has a ver3 5 4 tical ﬂat portion on the left. Fig. 5 shows this numerical range, as well as the numerical range of A1 for the unitary reduction of A to 1 A1 . Example 6. We will show that the following matrix has a numerical range whose boundary curve is of the same order, but the ﬂat portion is ‘‘hidden.’’

Fig. 5. Numerical ranges of the doubly stochastic matrix A and the primary component A1 in its reduction, having vertical ﬂat portions on the left.

46

K.A. Camenga et al. / Applied Mathematics and Computation 221 (2013) 40–47

0

0

B B1 B4 A¼B B5 B 12 @ 1 3

5 12

1 4

11 3

1 6

1 3

1C 4C

1 6

1 6

1C 4A

1 4

1 4

1 6

C C: C

The Hermitian part, H, again has three eigenspaces: the ﬁrst generated by e, with eigenvalue 1, the second generated by ½3; 1; 1; 1T with eigenvalue k2 ¼ 1 , and the third a 2-dimensional space generated again by ½0; 1; 0; 1T and ½0; 1; 1; 0T 3 3 and now corresponding to the double eigenvalue k3 ¼ 1 ¼ k4 . Computing H k3 I 1k eeT yields a symmetric negative semi 4 12 T deﬁnite matrix M of rank 1. Similarly to how this was done in Example 5, M ¼ ff , where f ¼ 4p1 ﬃﬃ3 ½3; 1; 1; 1T , and thus T 3 H ¼ k3 I þ 1k eeT ff . Furthermore, note that inequality (3.6) is satisﬁed: 1 1 11 . We conclude that when we reduce 4 12 12 12 Ato ð1Þ A1 , with

0

1 16 16 16 B C A1 ¼ @ 0 16 16 A; 1 0 0 6 @WðA1 Þ has a vertical ﬂat portion on the right. The convex hull of WðA1 Þ and f1g, which is the numerical range of A, does not have a vertical ﬂat portion. These numerical ranges are both illustrated in Fig. 6.

4. Numerical range of 6th degree All numerical ranges which do not fall into one of the above categories will have an ovular shaped boundary curve of degree 6. Example 7. In this example we will show that the following matrix does not satisfy the requirements of Theorems 1 and 3.

0 A¼

3 10 B3 B 10 B1 @ 10 3 10

1 10 1 5 2 5 3 10

3 10 3 10 1 5 1 5

3 1 10 1 C 5 C : 3 C A 10 1 5

We ﬁrst compute l ¼ 1=55 in (2.1). The characteristic polynomial of A evaluated at l is nonzero, and indeed the ﬁrst of the inequalities in (2.2) does not hold. Theorem 1 then implies that the numerical range of A is not the convex hull of an ellipse and two points. Next,

0 B 1 B H ¼ ðA þ AT Þ ¼ B @ 2

3 10 1 5 1 5 3 10

1 5 1 5 7 20 1 4

1 5 7 20 1 5 1 4

3 1 10 1 C 4 C C 1 A 4 1 5

pﬃﬃﬃﬃ 3 has four distinct eigenvalues (1, 20 , and 14017), which are therefore simple. Theorem 3 tells us that @WðAÞ does not contain arcs of a 4th degree curve. Altogether, we conclude based on our discussion at the end of Section 1 that WðAÞ is the convex hull of the point {1} and some ovular shape surrounded by a 6th degree boundary curve. Writing the reduction of A as ð1Þ A1 with

Fig. 6. The numerical range of A, and the numerical range of the component A1 which reveals a vertical ﬂat portion on the right.

K.A. Camenga et al. / Applied Mathematics and Computation 221 (2013) 40–47

47

Fig. 7. The numerical range of A, and the numerical range of the primary component in its reduction, A1 .

0

1 10

1 10 1 10

1 10

0

B A1 ¼ @ 0

1

1 5 1 C 10 A 1 10

tells us that this ovular shape actually is WðA1 Þ. Fig. 7 shows the numerical ranges of both A and A1 . References [1] [2] [3] [4] [5] [6] [7] [8] [9]

K.E. Gustafson, D.K.M. Rao, Numerical Range, The Field of Values of Linear Operators and Matrices, Springer, New York, 1997. R.A. Horn, C.R. Johnson, Topics in Matrix Analysis, Cambridge University Press, Cambridge, 1991. R. Kippenhahn, Über den Wertevorrat einer matrix, Math. Nachr. 6 (1951) 193–228. R. Kippenhahn, On the numerical range of a matrix, Linear Multilinear Algebra 56 (1–2) (2008) 185–225 (translated from the German by Paul F. Zachlin and Michiel E. Hochstenbach). D. Keeler, L. Rodman, I. Spitkovsky, The numerical range of 3 3 matrices, Linear Algebra Appl. 252 (1997) 115–139. L. Rodman, I.M. Spitkovsky, 3 3 matrices with a ﬂat portion on the boundary of the numerical range, Linear Algebra Appl. 397 (2005) 193–207. P. Rault, T. Sendova, I.M. Spitkovsky, 3-by-3 Matrices with elliptical numerical range revisited, Electron. Linear Algebra 26 (2013) 158–167. C.R. Johnson, An inclusion region for the ﬁeld of values of a doubly stochastic matrix based on its graph, Aequationes Math. 17 (2–3) (1978) 305–310. P. Nylen, T.Y. Tam, Numerical range of a doubly stochastic matrix, Linear Algebra Appl. 153 (1991) 161–176.