On the Distribution of a Scaled Condition Number

Alan Edelman* Department of Mathematics University of California Berkeley, California 94720

Abstract

In this note, we give the exact distribution of a scaled condition number used by Demmel to model the probability that matrix inversion is difficult. Specifically, consider a random matrix A and the scaled condition number κ_{D}(A) = ∥A∥_{F} · ∥A^{−1}∥. Demmel provided bounds for the condition number distribution when A has real or complex normally distributed elements. Here, we give the exact formula.

Key words and phrases: condition number, ill-conditioning, multivariate statistics, numerical analysis, random matrix ∗This work was started while the author was a Visiting Scholar in the Department of Mathematics, MIT, and completed at CERFACS with the support of the North Atlantic Treaty Organization under a Grant awarded in 1989.

1 Introduction and statement of results

In [4], Demmel investigates the probability that numerical analysis problems are difficult by unifying the common algebraic and geometric structures underlying the notion of ill-conditioning. As an application of his theory, he constructs a probabilistic model to examine the probability that matrix inversion is difficult. It is our goal in this note to work exclusively within his framework and derive exact distributions for the condition numbers that he considers. This condition number is a measure of difficulty in that the larger the value of the condition number, the more “difficult” matrix inversion becomes. The limitations of the model are discussed in [4]. We consider it a rather remarkable accident of mathematics that these distributions can be written down in a closed form at all. Although, as Demmel states, the assumption that matrices are uniformly distributed spherically is rather strong, the mathematics stands on its own, and indeed might have further applications to the “tubular neighborhoods” that he uses and perhaps also in multivariate statistics.

The objects of study are a scaled version κ_{D}(A) ≡ ∥A∥_{F} · ∥A^{−1}∥ of the usual condition number, and its distribution when considering random real and complex n by n matrices with elements distributed uniformly on the sphere 1 = ∥A∥_{F}^{2} = ∑ a_{ij}^{2}. Because of the scale invariance of the condition number and special properties of the normal distribution, it is equivalent to assume that the random matrices are generated with independent elements from a real or complex standard normal distribution. Demmel concludes that for real matrices

\frac{C(1 − 1/x)^{n^{2}−1}}{x} ≤ Prob(κ_{D}(A) ≥ x) ≤ ∑_{k=1}^{n^{2}} 2\binom{n^{2}}{k} (\frac{2n}{x})^{k} ,

where C is a constant. For complex matrices, on the other hand, he concludes that

\frac{(1 − x^{−1})^{2n^{2}−2}}{2n^{4}x^{2}} ≤ Prob(κ_{D}(A) ≥ x) ≤ \frac{e^{2}n^{5}(1 + n^{2}/x)^{2n^{2}−2}}{x^{2}} ,  (1.1)

and that asymptotically

Prob(κ_{D}(A) ≥ x) = \frac{n(n^{2} − 1)}{x^{2}} + o(\frac{1}{x^{2}} ),

as x → ∞ for fixed n.

In this note, we derive the exact probability distribution by combining exact distribution expressions for the smallest singular values [5, 6] of these random matrices with equations from [3, 8] which relate to κ_{D}. Our results are:

Probability Densities for κ_{D}(A): Real n by n matrices: μx^{1−n^{2}}(x^{2} − n)^{\frac{n(n+1)}{2} −2}_{2}F_{1} (\frac{n}{2} − \frac{1}{2}, \frac{n}{2} + 1; \frac{n^{2}}{2} + \frac{n}{2} − 1; −(x^{2} − n)) μ = \frac{2nΓ( \frac{n+1}{2} )Γ( \frac{n^{2}}{2} )}{\sqrt{π}Γ( \frac{n(n+1)}{2} −1)} Complex n by n matrices: 2n(n^{2} − 1)x^{1−2n^{2}}(x^{2} − n)^{n^{2}−2}

Here _{2}F_{1} denotes the Gauss hypergeometric function. Since obviously \sqrt{n} ≤ κ_{D}(A), these formulas are valid only for x ≥ \sqrt{n}.

For real matrices, the formula is cumbersome. For large x and n > 20, say, the formula below is quite adequate:

Prob(κ_{D}(A) ≥ x) ≈ n^{3/2}/x, x ≫ \sqrt{n}, n ≫ 1.

For complex matrices the exact distribution is a simple expression:

Prob(κ_{D}(A) ≥ x) = 1 − (1 − n/x^{2})^{n^{2}−1}, x > \sqrt{n}.

For large n, the condition numbers of real and complex matrices scale like n^{3/2}. To be precise let κ′ be the random variable 2κ_{D}/n^{3/2}. Then as n → ∞, for real matrices,

Prob(κ′ < x) → e^{−2/x−2/x^{2}}.  (1.2)

For complex matrices,

Prob(κ′ < x) → e^{−4/x^{2}}.  (1.3)

In fact, for large n the Demmel condition number κ_{D} of random uniformly distributed matrices is roughly \sqrt{n}/2 times as big as the ordinary 2-norm condition number κ_{2}. To be precise, as n → ∞, \frac{2}{\sqrt{n}}κ_{D}/κ_{2} converges almost surely to 1. Thus for large n, κ_{D} truly deserves to be called a scaled condition number, and the distribution of κ′ is the same as the distribution of κ_{2}/n which we have presented in [5, 6].

2 The distribution of κ_{D} (real case)

Let A be a real random n × n matrix with independent and identically distributed (iid) elements from a standard normal distribution. The matrix W = AA^{T} is said to be a Wishart matrix or to have the Wishart distribution. Our goal is to study the random quantity

κ_{D}(A) = \sqrt{\frac{∑_{i=1}^{n} λ_{i}}{λ_{n}}} ,

where λ_{1} ≥ … ≥ λ_{n} ≥ 0 are the eigenvalues of AA^{T}. Clearly κ_{D}(A) ≥ \sqrt{n}.

Let f_{n} be the probability density function (pdf) of (κ_{D}(A))^{−2} = λ_{n}/ ∑_{i=1}^{n} λ_{i}, and let g_{n} be the pdf of λ_{n}. The distribution function for κ_{D}(A) will be derived from two lemmas regarding f_{n} and g_{n}.

Lemma 2.1 (Davis) The pdfs f_{n} and g_{n} are related by

L ((1 + w)^{\frac{1}{2} n^{2}−2}f_{n} (\frac{1}{1 + w} )) (s) = 2Γ(n^{2}/2)e^{s}s^{− \frac{1}{2} n^{2}+1}g_{n}(2s),

where L denotes the Laplace transform.

This lemma was proved in [3], where the more general case of λ_{j}/ ∑_{i=1}^{n} λ_{i} is examined. These ratios arise in the multivariate analysis of variance (MANOVA) as described in multivariate analysis books such as [2].

The density function g_{n} is known exactly (see [5] or [6]):

Lemma 2.2 The density of the smallest eigenvalue of a Wishart matrix is

g_{n}(x) = \frac{n}{\sqrt{2π}} Γ (\frac{n + 1}{2}) x^{−1/2}e^{−xn/2}U(\frac{n − 1}{2}, −\frac{1}{2}, x/2).

When a > 0 and b < 1, the Tricomi function, U(a, b, z), is the unique solution to Kummer’s equation

z \frac{d^{2}w}{dz^{2}} + (b − z)\frac{dw}{dz} − aw = 0,  (2.1)

satisfying U(a, b, 0) = Γ(1 − b)/Γ(1 + a − b) and U(a, b, ∞) = 0.

Combining Lemma 2.1 and Lemma 2.2, we obtain

Theorem 2.1 The density of (κ_{D}(A))^{−2} is

f_{n}(x) = αx^{n^{2}/2−2}(x^{−1}−n)^{\frac{n(n+1)}{2} −2}_{2}F_{1} (\frac{n}{2} − \frac{1}{2}, \frac{n}{2} + 1; \frac{n^{2}}{2} + \frac{n}{2} − 1; −(x^{−1} − n)) ,

where

α = \frac{nΓ(\frac{n+1}{2})Γ(\frac{n^{2}}{2} )}{\sqrt{π}Γ(\frac{n(n+1)}{2} − 1)}

and _{2}F_{1} is the Gauss (hypergeometric) function.

Proof According to [7] (formula 7.522.4^{6}, p. 850),^{1} if b > 0,

L (w^{b−1}_{2}F_{1}(a, a − c + 1; b; −w)) (s) = Γ(b)s^{a−b}U(a, c, s),  (2.2)

where _{2}F_{1} is the Gauss (hypergeometric) function.

With a = \frac{n}{2} − \frac{1}{2}, b = \frac{n^{2}}{2} + \frac{n}{2} − 1, c = −\frac{1}{2}, and d = \frac{n}{2} + 1, we have from Lemmas 2.1 and 2.2 that

L ((1 + w)^{\frac{1}{2} n^{2}−2}f_{n} (\frac{1}{1 + w} )) (s) = αe^{s(1−n)}Γ(b)s^{a−b}U(a, c, s). (2.3)

From 2.2 we have

L (αw^{b−1}_{2}F_{1}(a, d; b; −w)) (s) = αΓ(b)s^{a−b}U(a, c, s).

Using familiar results concerning the Laplace transform, we then obtain

L (α(w − n + 1)^{b−1}_{2}F_{1}(a, d; b; −(w − n + 1))) (s) = αe^{s(1−n)}Γ(b)s^{a−b}U(a, c, s).

(2.4)

Combining 2.3 and 2.4, we obtain

(1 + w)^{\frac{1}{2} n^{2}−2}f_{n} (\frac{1}{1 + w} ) = α(w − n + 1)^{b−1}_{2}F_{1}(a, d; b; −(w − n + 1))

from which the theorem follows.

Corollary 2.1 Let h_{n} (x ≥ \sqrt{n}) be the density of κ_{D}(A) for real matrices. Then

h_{n}(x) = μx^{1−n^{2}}(x^{2} − n)^{\frac{n(n+1)}{2} −2}_{2}F_{1} (\frac{n}{2} − \frac{1}{2}, \frac{n}{2} + 1; \frac{n^{2}}{2} + \frac{n}{2} − 1; −(x^{2} − n)) ,

where

μ = \frac{2nΓ(\frac{n+1}{2})Γ(\frac{n^{2}}{2} )}{\sqrt{π}Γ(\frac{n(n+1)}{2} − 1)}

and _{2}F_{1} is the Gauss (hypergeometric) function.

1This formula is incorrect in older editions of [7]. We have verified that the formula as listed in our edition of [7] is indeed correct.

Proof This follows from Theorem 2.1 using the standard change of variable formula for probability densities.

Corollary 2.2 For fixed n, as x → ∞,

h_{n}(x) ∼ ν_{n}x^{−2},

where

ν_{n} = \frac{nΓ(\frac{n+1}{2})Γ(\frac{n^{2}}{2} )}{Γ(\frac{n^{2}−1}{2})Γ(\frac{n+2}{2})}.

For n > 20, ν_{n} ≈ n^{3/2}.

Proof The asymptotic formula for h_{n} follows from 15.3.4 and 15.1.20 of [1].

3 The distribution of κ_{D} (complex case)

The complex case is much easier than the real case, and the resulting formulas are considerably simpler. In [6], we gave a complete derivation of the exact density from first principles, but here we will proceed in an analogous manner to the real case.

Let A be a complex n × n matrix with independent and identically distributed (iid) elements from a complex standard normal distribution. A complex standard normal distribution can be defined as u + vi, where u and v are independent standard normals.

The matrix W = AA^{H} is said to be a complex Wishart matrix or have the complex Wishart distribution. Again our goal is to study the random quantity

κ_{D}(A) = \sqrt{\frac{∑_{i=1}^{n} λ_{i}}{λ_{n}}} ,

where λ_{1} ≥ … ≥ λ_{n} ≥ 0 are the eigenvalues of the complex Wishart matrix AA^{H}. As in the real case, κ_{D}(A) ≥ \sqrt{n}.

Using the same notation as in the real case, let f_{n} be the probability density function (pdf) of (κ_{D}(A))^{−2} = λ_{n}/ ∑_{i=1}^{n} λ_{i} and let g_{n} be the pdf of λ_{n}. The generalization of Lemma 2.1 for the complex case can be found in [8]. Lemma 3.1 (Krishnaiah and Schuurmann) The pdfs f_{n} and g_{n} are related by

L ((1 + w)^{n^{2}−2}f_{n} (\frac{1}{1 + w} )) (s) = Γ(n^{2})e^{s}s^{1−n^{2}}g_{n}(s).

Again, we have the density function g_{n} exactly (see [5] or [6]):

Lemma 3.2 The density of the smallest eigenvalue of a complex Wishart matrix is g_{n}(x) = ne^{−xn}, i.e, nλ_{min} is exponentially distributed.

Theorem 3.1 The density of (κ_{D}(A))^{−2} is

f_{n}(x) = n(n^{2} − 1)(1 − nx)^{n^{2}−2}.

Proof This formula can be derived from the two lemmas, and the integral formula for the gamma function.

Corollary 3.1 Let h_{n}(x) (x ≥ \sqrt{n}) be the density of the condition number κ_{D}(A) for complex matrices. Then

h_{n}(x) = 2n(n^{2} − 1)x^{1−2n^{2}}(x^{2} − n)^{n^{2}−2}.

Corollary 3.2 The probability distribution of κ_{D} is given in the complex case by

P (κ_{D} ≥ x) = 1 − (1 − n/x^{2})^{n^{2}−1}, x > \sqrt{n}.

The above result allows us to verify that indeed Corollary 3.3 For fixed n, as x → ∞,

P (κ_{D} ≥ x) ∼ n(n^{2} − 1)/x^{2}.

Acknowledgements

I would like to thank Jim Demmel for the inspiration for studying this condition number, and Tim Davis for the large amount of time he spent showing me how to use the Alliant FX/80 to perform random matrix experiments to check this work on eight processors in parallel. Also I wish to thank Nick Trefethen’s red pen for several useful suggestions.

References

[1] M. Abramowitz and I.A. Stegun, eds., Handbook of Mathematical Functions, Dover Publications, New York, 1970.

[2] T.W. Anderson, An Introduction to Multivariate Statistical Analysis, John Wiley & Sons, New York, 1958.

[3] A.W. Davis, On the ratios of the individual latent roots to the trace of a Wishart matrix, J. Mult. Anal. 2 (1972), 440–443.

[4] J.W. Demmel, The probability that a numerical analysis problem is difficult, Math. Comput. 50 (1988), 449–480.

[5] A. Edelman, Eigenvalues and condition numbers of random matrices, SIAM J. Matrix Anal. Appl. 9 (1988), 543–560.

[6] A. Edelman, Eigenvalues and Condition Numbers of Random Matrices, Ph.D. thesis, Dept. of Math., M.I.T., 1989.

[7] I.S. Gradshteyn and I.W. Ryzhik, Table of Integrals, Series, and Products, 6th printing, Academic Press, New York, 1980.

[8] P.R. Krishnaiah and F.J. Schuurmann, On the evaluation of some distributions that arise in simultaneous tests for the equality of the latent roots of the covariance matrix, J. Mult. Anal. 4 (1974), 265–282.

[9] F.J. Schuurmann, P.R. Krishnaiah, and A.K. Chattopadhyay, On the distributions of the ratios of the extreme roots to the trace of the Wishart matrix, J. Mult. Anal. 3 (1973), 445–453.