Matrices with Prescribed Row and Columns Sums-2

6. Random matrices with prescribed row and column sums Estimates of Theorems 2.1 and 3.1, however crude, allow us to obtain a description of a random or typical matrix from sets A0 (R, C) and A+ (R, C), considered as finite probability spaces with the uniform measures. Recall that x is a Bernoulli random variable if Pr {x = 0} = p and Pr {x = 1} = q for some p, q ≥ 0 such that p + q = 1. Clearly, E x = q. Recall that P0 (R, C) is the polytope of m × n matrices with row sums R, column sums C and entries between 0 and 1. Let function h : P0 (R, C) −→ R and the maximum entropy matrix Z0 ∈ P0 (R, C) be defined as in Lemma 5.2. The following result is proven in [Ba10a], see also [BH10a]. (6.1) Theorem. Suppose that polytope P0 (R, C) has a non-empty interior and let Z0 ∈ P0 (R, C) be the maximum entropy matrix. Let X = (xij ) be a random m × n matrix of independent Bernoulli random variables xij such that E X = Z0 . Then (1) The probability mass function of X is constant on the set A0 (R, C) of 0-1 matrices with row sums R and column sums C and Pr {X = D} = e−h(Z0 ) for all D ∈ A0 (R, C); (2) We have Pr {X ∈ A0 (R, C)} ≥ (mn)−γ(m+n) , where γ > 0 is an absolute constant. Theorem 6.1 implies that in many respects a random matrix D ∈ A0 (R, C) behaves as a random matrix X of independent Bernoulli random variables such that E X = Z0 , where Z0 is the maximum entropy matrix. More precisely, any event that is sufficiently rare for the random matrix X (that is, an event the probability of which is essentially smaller than (mn)−O(m+n) ), will also be a rare event for a random matrix D ∈ A0 (R, C). In particular, we can conclude that a typical matrix D ∈ A0 (R, C) is sufficiently close to Z0 as long as sums of entries over sufficiently large subsets S of indices are concerned. For an m × n matrix B = (bij ) and a subset n o S ⊂ (i, j) : i = 1, . . . , m, j = 1, . . . , n let σS (B) = X bij (i,j)∈S be the sum of the entries of B indexed by set S. We obtain the following corollary, see [Ba10a] for details. 20 (6.2) Corollary. Let us fix real numbers κ > 0 and 0 < δ < 1. Then there exists a number q = q(κ, δ) > 0 such that the following holds. Let (R, C) be margins such that n ≥ m > q and the polytope P0 (R, C) has a non-empty interior and let Z0 ∈ P0 (R, C) be the maximum entropy matrix. Let  S ⊂ (i, j) : i = 1, . . . , m; j = 1, . . . , n be a set such that σS (Z0 ) ≥ δmn and let ln ǫ = δ√ . m If ǫ ≤ 1 then n o Pr D ∈ A0 (R, C) : (1 − ǫ)σS (Z0 ) ≤ σS (D) ≤ (1 + ǫ)σS (Z0 ) ≥ 1 − n−κn . Recall that x is a geometric random variable if Pr {x = k} = pq k for k = 0, 1, 2, . . . for some p, q ≥ 0 such that p + q = 1. We have E x = q/p. Recall that P+ (R, C) is the polytope of m × n non-negative matrices with row sums R and column sums C. Let function g : P+ (R, C) −→ R and the maximum entropy matrix Z+ ∈ P0 (R, C) be defined as in Lemma 5.3. The following result is proven in [Ba10b], see also [BH10a]. (6.3) Theorem. Let Z+ ∈ P0 (R, C) be the maximum entropy matrix. Let X = (xij ) be a random m×n matrix of independent geometric random variables xij such that E X = Z+ . Then (1) The probability mass function of X is constant on the set A+ (R, C) of nonnegative integer matrices with row sums R and column sums C and Pr {X = D} = e−g(Z+ ) for all D ∈ A+ (R, C); (2) We have Pr {X ∈ A+ (R, C)} ≥ N −γ(m+n) , where γ > 0 is an absolute constant and N = r1 + . . . + rm = c1 + . . . + cn for R = (r1 , . . . , rm ) and C = (c1 , . . . , cn ). Theorem 6.3 implies that in many respects a random matrix D ∈ A+ (R, C) behaves as a matrix X of independent geometric random variables such that E X = Z+ , where Z+ is the maximum entropy matrix. More precisely, any event that is sufficiently rare for the random matrix X (that is, an event the probability of which is essentially smaller than N −O(m+n) ), will also be a rare event for a random matrix D ∈ A+ (R, C). In particular, we can conclude that a typical matrix D ∈ A+ (R, C) is sufficiently close to Z+ as long as sums of entries over sufficiently large subsets S of indices are concerned. Recall that σS (B) denotes the sum of the entries of a matrix B indexed by a set S. We obtain the following corollary, see [Ba10b] for details. 21 (6.4) Corollary. Let us fix real numbers κ > 0 and 0 < δ < 1. Then there exists a positive integer q = q(κ, δ) such that the following holds. Let R = (r1 , . . . , rm ) and C = (c1 , . . . , cm ) be positive integer vectors such that r1 + . . . + rm = c1 + . . . + cn = N , N δN ≤ ri ≤ for i = 1, . . . m, m δm δN N ≤ cj ≤ for j = 1, . . . , n n δn and N ≥ δ. mn  Suppose that n ≥ m > q and let S ⊂ (i, j) : i = 1, . . . , m, j = 1, . . . , n be a set such that |S| ≥ δmn. Let Z+ ∈ P+ (R, C) be the maximum entropy matrix and let ln n ǫ = δ 1/3 . m If ǫ ≤ 1 then n o Pr D ∈ A+ (R, C) : (1 − ǫ)σS (Z+ ) ≤ σS (D) ≤ (1 + ǫ)σS (Z+ ) ≥ 1 − n−κn . As is discussed in [BH10a], the ultimate reason why Theorems 6.1 and 6.3 hold true is that the matrix X of independent Bernoulli random variables such that E X = Z0 is the random matrix with the maximum possible entropy among all random m × n matrices with 0-1 entries and the expectation in the affine subspace of the matrices with row sums R and column sums C and the matrix X of independent geometric random variables such that E X = Z+ is the random matrix with the maximum possible entropy among all random m × n matrices with non-negative integer entries and the expectation in the affine subspace of the matrices with row sums R and column sums C. Thus Theorems 6.1 and 6.3 can be considered as an illustration of the Good’s thesis [Go63] that the “null hypothesis” for an unknown probability distribution from a given class should be the hypothesis that the unknown distribution is, in fact, the distribution of the maximum entropy in the given class. (6.5) Sketch of proof of Theorem 6.1. Let Z0 = (zij ) be the maximum entropy matrix as in Lemma 5.2. Let us choose D ∈ A0 (R, C), D = (dij ). Using (5.2.1), 22 we get Y d Y e(λi +µj )dij  1−d Pr X = D = zijij (1 − zij ) ij = 1 + eλi +µj i,j ij   m n X Y X 1 = exp λi ri + µj cj   1 + eλi +µj i=1 j=1 ij =e−h(Z0 ) , which proves Part (1). To prove Part (2), we use Part (1), Theorem 2.1 and Lemma 5.2. We have  Pr X ∈ A0 (R, C) = |A0 (R, C)| e−h(Z0 ) ≥ (mn)−γ(m+n) α0 (R, C)e−h(Z0 ) =(mn)−γ(m+n) for some absolute constant γ > 0.  (6.6) Sketch of proof of Theorem 6.3. Let Z+ = (zij ) be the maximum entropy matrix as in Lemma 5.3. Let us choose D ∈ A+ (R, C), D = (dij ). Using (5.3.1), we get   dij Y
1 zij 1 − e−λi −µj e−(λi +µj )dij = 1 + zij 1 + zij i,j ij   n m Y  X X
µj cj λi ri − = exp − 1 − e−λi −µj   Y  Pr X = D = i=1 j=1 ij =e−g(Z+ ) , which proves Part (1). To prove Part (2), we use Part (1), Theorem 3.1 and Lemma 5.3. We have  Pr X ∈ A+ (R, C) = |A+ (R, C)| e−g(Z+ ) ≥ N −γ(m+n) α+ (R, C)e−g(Z+ ) =N −γ(m+n) for some absolute constant γ > 0.  (6.7) Open questions. Theorems 6.1 and 6.3 show that a random matrix D ∈ A0 (R, C), respectively D ∈ A+ (R, C), in many respects behaves like a matrix of independent Bernoulli, respectively geometric, random variables whose expectation is the maximum entropy matrix Z0 , respectively Z+ . One can ask whether individual entries dij of D behave asymptotically as Bernoulli, respectively geometric, random variables with expectations zij as the size of the matrices grows. In the simplest situation we ask the following 23 (6.7.1) Question. Let (R, C) be margins and let (Rk , Ck ) be margins obtained from (R, C) by cloning as in Section 5.5. Is it true that as k grows, the entry d11 of a random matrix D ∈ A0 (Rk , Ck ), respectively D ∈ A+ (Rk , Ck ), converges in distribution to the Bernoulli, respectively geometric, random variable with expectation z11 , where Z0 = (zij ), respectively Z+ = (zij ), is the maximum entropy matrix of margins (R, C)? Some entries of the maximum entropy matrix Z+ may turn out to be surprisingly large, even for reasonably looking margins. In [Ba10b], the following example is considered. Suppose that m = n and let Rn = Cn = (3n, n, . . . , n). It turns out that the entry z11 of the maximum entropy matrix Z+ is linear in n, namely z11 > 0.58n, while all other entries remain bounded by a constant. One can ask whether the d11 entry of a random matrix D ∈ A+ (Rn , Cn ) is indeed large, as the value of z11 suggests. (6.7.2) Question. Let (Rn , Cn ) be margins as above. Is it true that as n grows, one has E d11 = Ω(n) for a random matrix D ∈ A+ (Rn , Cn )? Curiously, the entry z11 becomes bounded by a constant if 3n is replaced by 2n. 7. Asymptotic formulas for the number of matrices with prescribed row and column sums In this section, we discuss asymptotically exact estimates for |A0 (R, C)| and |A+ (R, C)|. (7.1) An asymptotic formula for |A0 (R, C)|. Theorem 6.1 suggests the following way to estimate the number |A0 (R, C)| of 0-1 matrices with row sums R and column sums C. Let us consider the matrix of independent Bernoulli random variables as in Theorem 6.1 and let Y be the random (m + n)-vector obtained by computing the row and column sums of X. Then, by Theorem 6.1, we have   (7.1.1) |A0 (R, C)| = eh(Z0 ) Pr X ∈ A0 (R, C) = eh(Z0 ) Pr Y = (R, C) . Now, random (m + n)-vector Y is obtained as a sum of mn independent random  vectors and E Y = (R, C), so it is not unreasonable to assume that Pr Y = (R, C) can be estimated via some version of the Local Central Limit Theorem. In [BH10b] we show that this is indeed the case provided one employs the Edgeworth correction factor in the Central Limit Theorem. We introduce the necessary objects to state the asymptotic formula for the number of 0-1 matrices with row sums R and column sums C. Let Z0 = (zij ) be the maximum entropy matrix as in Lemma 5.2. We assume that 0 < zij < 1 for all i and j. Let us consider the quadratic form q0 : Rm+n −→ R defined by
1 X 2 2 zij − zij (si + tj ) q0 (s, t) = 2 1≤i≤m 1≤j≤n for s = (s1 , . . . , sm ) 24 and t = (t1 , . . . , tn ) . Quadratic form q0 is positive semidefinite with the kernel spanned by vector   u = 1, . . . , 1; −1, . . . , −1 . {z } | {z } | m times n times Let H = u⊥ be the hyperplane in Rm+n defined by the equation (7.1.2) s 1 + . . . + s m = t1 + . . . + tn . Then the restriction q0 |H of q0 onto H is a positive definite quadratic form and we define its determinant det q0 |H as the product of the non-zero eigenvalues of q0 . We consider the Gaussian probability measure on H with the density proportional to e−q0 and define random variables φ0 , ψ0 : H −→ R by φ0 (s, t) = 1 X 3 zij (1 − zij ) (2zij − 1) (si + tj ) 6 and 1≤i≤m 1≤j≤n ψ0 (s, t) =
1 X 4 2 zij (1 − zij ) 6zij − 6zij + 1 (si + tj ) 24 1≤i≤m 1≤j≤n for (s, t) = (s1 , . . . , sm ; t1 , . . . , tn ) . We let µ0 = E φ20 and ν0 = E ψ0 . (7.2) Theorem. Let us fix 0 < δ < 1/2, let R = (r1 , . . . , rm ) and C = (c1 , . . . , cn ) be margins such that m ≥ δn and n ≥ δm. Let Z0 = (zij ) be the maximum entropy matrix as in Lemma 5.2 and suppose that δ ≤ zij ≤ 1 − δ for all i and j. Let the quadratic form q0 and values µ0 and ν0 be as defined in Section 7.1. Then the number √ n µ o eh(Z0 ) m + n 0 p (7.2.1) exp − + ν0 2 (4π)(m+n−1)/2 det q0 |H approximates the number |A0 (R, C)| of as m, n −→ +∞ within a relative error which approaches 0 as m, n −→ +∞. More precisely, for any 0 < ǫ ≤ 1/2, the value of (7.2.1) approximates |A0 (R, C)| within relative error ǫ provided  γ(δ) 1 m, n ≥ ǫ for some γ(δ) > 0. Some remarks are in order. 25 All the ingredients of formula (7.2.1) are efficiently computable, in time polynomial in m + n, see [BH10b] for details. If all row sums are equal then we have zij = cj /m by symmetry and if all column sums are equal, we have zij = ri /n. In particular, if all row sums are equal and if all column sums are equal, we obtain the asymptotic formula of [C+08]. Let us consider formula (7.1.1). If, in the spirit  of the Local Central Limit  Theorem, we approximated Pr Y = (R, C) by Pr Y ∗ ∈ (R, C) + Π , where Y ∗ is the (m + n − 1)-dimensional random Gaussian vector whose expectation and covariance matrix match those of Y and where Π is the set of points on the hyperplane H that are closer to (R, C) than to any other integer vector in H, we would have obtained the first part √ eh(Z0 ) m + n p (4π)(m+n−1)/2 det q0 |H of formula (7.2.1). Under the conditions of Theorem 7.2 we have n µ o 0 c1 (δ) ≤ exp − + ν0 ≤ c2 (δ) 2 for some constants c1 (δ), c2 (δ) > 0 and this factor represents the Edgeworth correction to the Central Limit Theorem. We note that the constraints δ ≤ zij ≤ 1−δ are, generally speaking, unavoidable. If the entries zij of the maximum entropy matrix are uniformly small, then the distribution of the random vector Y of row and column sums of the random Bernoulli matrix X is no longer approximately Gaussian but approximately Poisson and formula (7.2.1) does not give correct asymptotics. The sparse case of small row and column sums is investigated in [G+06]. More generally, to have some analytic formula approximating |A0 (R, C)| we need certain regularity conditions on (R, C), since the number |A0 (R, C)| becomes volatile when the margins (R, C) approach the boundary of the Gale-Ryser conditions, cf. [JSM92]. By requiring that the entries of maximum entropy matrix Z0 are separated from both 0 and 1, we ensure that the margins (R, C) remain sufficiently inside the polyhedron defined by the Gale-Ryser inequality and the number of 0-1 matrices with row sums R and column sums C changes sufficiently smoothly when R and C change. (7.3) An asymptotic formula for |A+ (R, C)|. As in Theorem 6.3, let X be the matrix of independent geometric random variables such that E X = Z+ , where Z+ is the maximum entropy matrix. Let Y be the random (m + n)-vector obtained by computing the row and column sums of X. Then, by Theorem 6.3, we have (7.3.1)   |A+ (R, C)| = eg(Z+ ) Pr X ∈ A+ (R, C) = eg(Z+ ) Pr Y = (R, C) . In [BH09] we show how to estimate the probability that Y = (R, C) using the Local Central Limit Theorem with the Edgeworth correction. 26 Let Z+ = (zij ) be the maximum entropy matrix as in Lemma 5.3. Let us consider the quadratic form q+ : Rm+n −→ R defined by q+ (s, t) = 1 X 2 1≤i≤m 1≤j≤n for
2 2 zij + zij (si + tj ) s = (s1 , . . . , sm ) and t = (t1 , . . . , tn ) . Let H ⊂ Rm+n be the hyperplane defined by (7.1.2). The restriction q+ |H of q+ onto H is a positive definite quadratic form and we define its determinant det q+ |H as the product of the non-zero eigenvalues of q+ . We consider the Gaussian probability measure on H with the density proportional to e−q+ and define random variables φ+ , ψ+ : H −→ R by φ+ (s, t) = 1 X 3 zij (1 + zij ) (2zij + 1) (si + tj ) 6 and 1≤i≤m 1≤j≤n ψ+ (s, t) =
1 X 4 2 zij (1 + zij ) 6zij + 6zij + 1 (si + tj ) 24 1≤i≤m 1≤j≤n for (s, t) = (s1 , . . . , sm ; t1 , . . . , tn ) . We let µ+ = E φ2+ and ν+ = E ψ+ . (7.4) Theorem. Let us fix 0 < δ < 1, let R = (r1 , . . . , rm ) and C = (c1 , . . . , cn ) be margins such that m ≥ δn and n ≥ δm. Let Z+ = (zij ) be the maximum entropy matrix as in Lemma 5.3. Suppose that δτ ≤ zij ≤ τ for all i, j for some τ ≥ δ. Let the quadratic form q+ and values µ+ and ν+ be as defined in Section 7.3. Then the number √ o n µ eg(Z+ ) m + n + p + ν+ (7.4.1) exp − 2 (4π)(m+n−1)/2 det q+ |H approximates the number |A+ (R, C)| of as m, n −→ +∞ within a relative error which approaches 0 as m, n −→ +∞. More precisely, for any 0 < ǫ ≤ 1/2, the value of (7.4.1) approximates |A+ (R, C)| within relative error ǫ provided  γ(δ) 1 m, n ≥ ǫ 27 for some γ(δ) > 0. All the ingredients of formula (7.4.1) are efficiently computable, in time polynomial in m + n, see [BH09] for details. If all row sums are equal then we have zij = cj /m by symmetry and if all column sums are equal, we have zij = ri /n. In particular, if all row sums are equal and if all column sums are equal, we obtain the asymptotic formula of [C+07]. The term √ eg(Z+ ) m + n p (4π)(m+n−1)/2 det q+ |H corresponds to the Gaussian approximation for the distribution of the random vector Y in (7.3.1), while n µ o + exp − + ν+ 2 is the Edgeworth correction factor. While the requirement that the entries of the maximum entropy matrix Z+ are separated from 0 is unavoidable (if zij are small, the coordinates of Y are asymptotically Poisson, not Gaussian, see [GM08] for the analysis of the sparse case), it is not clear whether the requirement that all zij are within a constant factor of each other is indeed needed. It could be that around certain margins (R, C) the number |A+ (R, C)| experiences sudden jumps, as the margins change, which precludes the existence of an analytic expression similar to (7.4.1) for |A+ (R, C)|. A candidate for such an abnormal behavior is supplied by the margins discussed in Section 6.7. Namely, if m = n and R = C = (λn, n, . . . , n) then for λ = 2 all the entries of the maximum entropy matrix Z+ are O(1), while for λ = 3 the first entry z11 grows linearly in n. Hence for some particular λ between 2 and 3 a certain “phase transition”’ occurs: the entry z11 jumps from O(1) to Ω(n). It would be interesting to find out if there is indeed a sharp change in |A+ (R, C)| when λ changes from 2 to 3. 8. Concluding remarks Method of Sections 6 and 7 have been applied to some related problems, such as counting higher-order “tensors” with 0-1 or non-negative integer entries and prescribed sums along coordinate hyperplanes [BH10a] and counting graphs with prescribed degrees of vertices [BH10b], which corresponds to counting symmetric 0-1 matrices with zero trace and prescribed row (column) sums. In general, the problem can be described as follows: we have a polytope P ⊂ Rd defined as the intersection of the non-negative orthant Rd+ with an affine subspace A in Rd and we construct a d-vector X of independent Bernoulli (in the 0-1 case) or geometric (in the non-negative integer case) random variables, so that the expectation of X lies in A and the distribution of X is uniform, when restricted onto the set of 0-1 or integer points in P . Random vector X is determined by its expectation E X = z and z is found by solving a convex optimization problem on 28 P . Since vector X conditioned on the set of 0-1 or non-negative integer vectors in P is uniform, the number of 0-1 or non-negative integer points in P is expressed in terms of the probability that X lies in A. Assuming that the affine subspace A is defined by a system Ax = b of linear equations, where A is k × d matrix of rank k < d, we define a k-vector Y = AX of random variables and estimate the probability that Y = b by using a Local Central Limit Theorem type argument. Here we essentially use that E Y = b, since the expectation of X lies in A. Not surprisingly, the argument works the easiest when the codimension k of the affine subspace (and hence the dimension of vector Y ) is small. In particular, counting higher-order “tensors” is easier than counting matrices, the need in the Edgeworth correction factor, for example, disappears as the vector Y turns out to be closer in distribution to a Gaussian vector, see  [BH10a]. Once a Gaussian or almost Gaussian estimate for the probability Pr Y = b is established, one can claim a certain concentration of a random 0-1 or integer point in P around z = E X. References [Ba07] [Ba09] [Ba10a] [Ba10b] [BH09] [BH10a] [BH10b] [Be74] [BR91] [C+08] [C+07] [Ga02] [Go63] [Go76] A. Barvinok, Brunn-Minkowski inequalities for contingency tables and integer flows, Advances in Mathematics 211 (2007), 105–122. A. Barvinok, Asymptotic estimates for the number of contingency tables, integer flows, and volumes of transportation polytopes, International Mathematics Research Notices 2009 (2009), 348–385. A. Barvinok, On the number of matrices and a random matrix with prescribed row and column sums and 0-1 entries, Advances in Mathematics 224 (2010), 316–339. A. Barvinok, What does a random contingency table look like?, Combinatorics, Probability and Computing 19 (2010), 517–539. A. Barvinok and J.A. Hartigan, An asymptotic formula for the number of non-negative integer matrices with prescribed row and column sum, preprint arXiv:0910.2477 (2009). A. Barvinok and J.A. Hartigan, Maximum entropy Gaussian approximation for the number of integer points and volumes of polytopes, Advances in Applied Mathematics 45 (2010), 252–289. A. Barvinok and J.A. Hartigan, The number of graphs and a random graph with a given degree sequence, preprint arXiv:1003.0356 (2010). E. Bender, The asymptotic number of non-negative integer matrices with given row and column sums, Discrete Math. 10 (1974), 217–223. R.A. Brualdi and H.J. Ryser, Combinatorial Matrix Theory, Encyclopedia of Mathematics and its Applications, 39, Cambridge University Press, Cambridge, 1991. E.R. Canfield, C. Greenhill, and B.D. McKay, Asymptotic enumeration of dense 01 matrices with specified line sums, Journal of Combinatorial Theory. Series A 115 (2008), 32–66. E.R. Canfield and B.D. McKay, Asymptotic enumeration of contingency tables with constant margins, preprint arXiv math.CO/0703600, Combinatorica, to appear (2007). R.J. Gardner, The Brunn-Minkowski inequality, Bull. Amer. Math. Soc. (N.S.) 39 (2002), 355–405. I.J. Good, Maximum entropy for hypothesis formulation, especially for multidimensional contingency tables, Ann. Math. Statist. 34 (1963), 911–934. I.J. Good, On the application of symmetric Dirichlet distributions and their mixtures to contingency tables, Ann. Statist. 4 (1976), 1159–1189. 29 [GC77] [G+06] [GM08] [JSM92] [Kh57] [LW01] [NN94] I.J. Good and J.F. Crook, The enumeration of arrays and a generalization related to contingency tables, Discrete Mathematics 19 (1977), 23–45. C. Greenhill, B.D. McKay, and X. Wang, Asymptotic enumeration of sparse 0-1 matrices with irregular row and column sums, Journal of Combinatorial Theory. Series A 113 (2006), 291–324. C. Greenhill and B.D. McKay, Asymptotic enumeration of sparse nonnegative integer matrices with specified row and column sums, Advances in Applied Mathematics 41 (2008), 59–481. M. Jerrum, A. Sinclair and B. McKay, When is a graphical sequence stable?, Random Graphs, Vol. 2 (Poznań, 1989), Wiley-Intersci. Publ., Wiley, New York, 1992, pp. 101– 115. A.I. Khinchin, Mathematical Foundations of Information Theory, Dover Publications Inc., New York, N. Y., 1957. J.H. van Lint and R.M. Wilson, A Course in Combinatorics. Second edition, Cambridge University Press, Cambridge, 2001. Y. Nesterov and A. Nemirovskii, Interior-Point Polynomial Algorithms in Convex Programming, SIAM Studies in Applied Mathematics, 13, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1994. 30