Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 94
Текст из файла (страница 94)
See Table 26.1 for the first few conditionnumbers (these were obtained by computing the inverse exactly using MATLAB’S Symbolic Math Toolbox [204, 1993] and then computing the norm ofthe numeric representation of the inverse; the numbers given are correct tothe figures shown).It is an interesting fact that the matrix= ( 1 / (i + j) )(asubmatrix of H n+l ) satisfies µ n :== π + O(1/log n) asasproved by Taussky [996, 1949]. The convergence to π is very slow: µ 200 = 2.01,µ300 = 2.08, µ400 = 2.12.That H n -1 is known explicitly is not as useful a property for testing aninversion algorithm as it might appear, because Hn cannot be stored exactlyin floating point arithmetic.
This means that when we attempt to invert H nwe at best invert H n + ∆H (the matrix actually stored on the computer),A G ALLERY516OFT EST M ATRICESTable 26.1. Condition numbers of Hilbert and Pascal matrices.n2345678910111213141516K2(Hn)1.9281e15.2406e21.5514e44.7661e51.4951e74.7537e81.5258e104.9315e111.6026e135.2307e141.7132e165.6279e171.8534e196.l166e202.0223e22K2(Pn)6.8541e06.1984el6.9194e28.5175e31.1079e51.4934e62.0645e72.9078e84.1552e96.0064e108.7639e111.2888e131.9076e142.8396e154.2476e16w h e r e | ∆H| < uH n , and (Hn + ∆H)–1 can differ greatly from Hn-1, in viewof the ill conditioning.
A possible way round this difficulty is to start withthe integer matrix H n -1 , but its entries are so large that they are exactlyrepresentable in IEEE double precision arithmetic only for n less than 13.The Hilbert matrix is a special case of a Cauchy matrix C nwhose elements are cij = 1/(xi + yj), where x, yare given n-vectors(take xi = yi = i – 0.5 for the Hilbert matrix). The following formulae give theinverse and determinant of Cn, and therefore generalize those for the Hilbertmatrix. The (i, j) element of Cn-1 isandthe latter formula having been published by Cauchy in 1841 [189, 1841,pp.
151–159]. The LDU factors of Cn have been found explicitly by Gohberg26.2 R ANDOM M ATRICES517and Koltracht [454, 1990]: lkk = ukk = 1 andIt is known that Cn is totally positive if 0 < x1 < · · · < xn and 0 < yl <· · · < yn (as is true for the Hilbert matrix) [998, 1962, p. 295]. Interestingly,the sum of all the elements of[667, 1973, Ex. 44, §1.2.3].26.2. Random MatricesRandom matrices are widely used for test purposes.
Perhaps the earliest useof random matrices in numerical analysis was by von Neumann and Goldstine [1057, 1947], [462, 1951], who estimated the size of their error bounds formatrix inversion (see §9.6) for the case of random matrices; to do this, theyhad to estimate the condition number of a random matrix.Intuitively, one might expect random matrices to be good at revealingprogramming errors and unusual behaviour of algorithms, but this expectation is not necessarily correct.
For example, Miller [759, 1984, pp. 96-97]describes a mutation experiment involving Fortran codes for Gaussian elimination without pivoting, Gaussian elimination with partial pivoting, and Gauss–Jordan elimination with partial pivoting. For each code, all possible mutantswere generated, where a mutant is obtained by making a single typographicalchange to the source code. All the mutants were tested on a single randomlinear system Ax = b, with known solution, where a ij was chosen from theuniform [– 1, 1] distribution. Many mutants were detected by their failure topass the test of producing a solution with forward error less than a tolerance.However, some mutants passed this test, including all those that solve a system correctly in exact arithmetic; mutants in the latter class include thosethat select an incorrect pivot row and thus implement a numerically unstablealgorithm.
A conclusion to be drawn from Miller’s experiment is that randomtest data can reveal some programming errors, but will not reveal all.A good example of a problem for which random test matrices are very poorat revealing algorithmic weaknesses is matrix condition number estimation.The popular condition estimation algorithms can yield poor estimates but, in518A G ALLERYOFT EST M ATRICESpractice, never produce them for a random matrix (see Chapter 14). The roleof random matrices here is to indicate the average quality of the estimates.Edelman [340, 1993] summarizes the properties of random matrices wellwhen he says thatWhat is a mistake is to psychologically link a random matrix withthe intuitive notion of a “typical” matrix or the vague concept of“any old matrix.” In contrast, we argue that “random matrices”are very special matrices.
The larger the size of the matrices themore predictable they are because of the central limit theorem.Various results are known about the behaviour of matrices with elementsfrom the normal N(0, 1) distribution. Matrices of this type are generated byMATLAB’S randn function. Let An denote an n x n matrix from this distribution and let E(·) be the expectation operator. Then, in the appropriateprobabilistic sense, the following results hold as n(real data),(26.5)(complex data),(26.6)(real data),(complex data),(real or complex data).(26.7)(26.8)(26.9)For details of (26.5)-(26.8) see Edelman [335, 1988].
Edelman conjectures thatthe condition number results are true for any distribution with mean O—inparticular, the uniform [– 1, 1] distribution used by MATLAB’S rand function.The results (26.5) and (26.6) show that random matrices from the normalN(0, 1) distribution tend to be very well conditioned.The spectral radius result (26.9) has been proved as an inequality by Geman [432, 1986] for independent identically distributed random variables aijwith zero mean and unit variance, and computer experiments suggest theapproximate equality for the standard normal distribution [432, 1986].A question of interest in eigenvalue applications is how many eigenvaluesof a random real matrix are real. The answer has been given by Edelman,Kostlan, and Shub [344, 1994]: denoting by En the expected number of realeigenvalues of an n x n matrix from the normal N(0, 1) distribution,Thus the proportion of real eigenvalues, En/n, tends to zero as nformulae for En for finite n are also given in [344, 1994].Exact26.3 “R ANDSVD ” M ATRICES51926.3.
“Randsvd” MatricesBy randsvd25 we mean a matrix Aformed as A = UΣVT, whereare random orthogonal matrices and Σ = diag(σ i )is a given matrix of singular values. This type of matrix has apredetermined singular value distribution (and 2-norm condition number),but is, nevertheless, random.Randsvd matrices have been widely used as test matrices, for examplefor condition estimators [534, 1987], [946, 1980], and in the LAPACK testingsoftware. Singular value distributions of interest include, with α := K2(A) > 1a parameter,1. one large singular value: σ1 = 1, σi = α–1, i = 2:n;2.
one small singular value: σi = 1, i = 1:n – 1, σ n = α– l ;3. geometrically distributed singular values: σi = β 1–i,i=1:n, where4. arithmetically distributed singular values: σ i = 1 – (1 – α - l )(i – 1)/(n –1), i = 1:n.To be precise about what we mean by “random orthogonal matrix” wespecify matrices from the Haar distribution, which is a natural distributionover the space of orthogonal matrices, defined in terms of a measure calledthe Haar measure [780, 1982, 52.1.4].
If Ahas elements from thenormal N(0, σ 2) distribution and A has the QR factorization A = QR, withthe factorization normalized so that the diagonal elements of R are positive,then Q is from the Haar distribution, for any variance σ 2 [100, 1979], [946,1980]. If we compute Q from this prescription, the cost is 2n3 flops. For ourrandsvd application, a more efficient approach is based on the following resultof Stewart [946, 1980].Theorem 26.1 (Stewart). Let the independent vectors x ihave elements from the normal N(0, 1) distribution for i = 1:n – 1. Let Pi =diagwhereis the Householder transformation that reduces xito rii el. Then the product Q = DP1P2 . . . Pn-l is a random orthogonal matrix from the Haar distribution, where D = diag(sign(rii ) ) .This result allows us to compute a product form representation of a randomn x n orthogonal matrix from the Haar distribution in O(n2 ) flops.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
