Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 95
Текст из файла (страница 95)
If weimplicitly compute Uandusing the construction inTheorem 26.1, and then form A = UΣVT, exploiting the structure, the cost25Randsvd is the name of the MATLAB M-file in the Test Matrix Toolbox that generatesmatrices of this type.A G ALLERY520OFT EST M ATRICESis m3 + n3 flops. The alternative, which involves obtaining U and V from theQR factorization of random matrices A, is about twice as expensive.Note that forming UΣVT with U and V single Householder matrices, as issometimes done in the literature, is not recommended, as it produces matricesof a very special form: diagonal plus a rank-2 correction.The construction of randsvd matrices is, of course, easily adapted to produce random symmetric matrices A = QΛQT with given eigenvalues.26.4.
The Pascal MatrixThe numbers in Pascal’s triangle satisfy, practically speaking,infinitely many identities.— RONALD L. GRAHAM, DONALD E. KNUTH, and OREN PATASHNIK,Concrete Mathematics (1989)A particularly entertaining test matrix is the Pascal matrix Pndefined byThe rows of Pascal’s triangle appear as anti-diagonals of Pn26:>> p = pascal(6)P=111111123456136101521141020355615153570126162156126252The earliest references to the Pascal matrix appear to be in 1958 by Newmanand Todd [796, 1958] and by Rutishauser [887, 1958] (see also Newman [795,1962, pp. 240-24 1]); Newman and Todd say that the matrix was introducedto them by Rutishauser. The matrix was independently suggested as a testmatrix by Caffney [177, 1963].Rutishauser [886, 1968, §8] notes that Pn belongs to the class of momentmatrices M whose elements are contour integrals26In the MATLAB displays below we use the pascal function from the Test Matrix Toolbox.This differs from the pascal function supplied with Matlab 4.2 only in that pascal (n, 2) isrearranged.52126.4 T HE P ASCAL M ATRIXAll moment matrices corresponding to a positive weight function w(z) on acontour C are Hermitian positive definite (as is easily verified by consideringthe quadratic form y*My).
The choice C = [0, 1], with weight function w(z) =1, yields the Hilbert matrix. The Pascal matrix is obtained for C the circle{z : |z – 1| = 1 } and w(z) = (2 πi(z – 1))-1 (not w(z) = (2π)-1 as stated in[886, 1968]); the change of variable z = 1 + exp(iθ ) yields a moment integralwith a positive weight function.Remarkably, the Cholesky factor of the Pascal matrix again contains therows of Pascal’s triangle, now arranged columnwise:>> R = chol(P)R=10000011000012100013310014641015101051The scaled and transposed Cholesky factor S = RT diag(1, –1,1, –1,. . . . (–1)n + 1 )is returned by pascal (n, 1 ):>> S = pascal(6, 1)S =1111110-1-2-3-4-50013610000-1-4-1000001500000-1It is involuntary: S 2 = I.
This special property leads us to several moreproperties of P = P n . First, since P = SS T , P –1 = S -T S -l = S T S, and soP-1 has integer entries (as is also clear from the fact that det(P) = det(R)2 =1). Moreover,P = SST = S(STS)S-1 = SP-1S-1,so P and P–1 are similar and hence have the same eigenvalues. In otherwords, the eigenvalues appear in reciprocal pairs. In fact, the characteristicpolynomial π n has a palindromic coefficient vector, which implies the reciprocal eigenvalues property, since π n (λ) =This is illustrated asfollows (making use of MATLAB’S Symbolic Math Toolbox):522A G ALLERYOFT EST M ATRICES>> charpoly(P)ans =1-351*x+6084*x^2-13869*x^3+6084*x^4-351*x^5+x^6>> eig(P)ans =0.00300.06430.48932.043615.5535332.8463Since P is symmetric, its eigenvalues are its singular values and so we alsohave that ||P||2 = ||P-1 ||2 and ||P||F = ||P-1 ||F.
Nowwhere for the last equality we used a binomial coefficient summation identity from [477, 1989, p. 161]. Hence, using Stirling’s approximation (n ! ~Thus Pn is exponentially ill conditioned as nIt is worth pointing out that it is not hard to generate symmetric positivedefinite matrices with determinant 1 and the reciprocal root property. LetX = ZDZ -l where Z is nonsingular and D = diag(±1) ±I. Then X2 = Iand the matrix A = X T X has the desired properties. If we choose Z lowertriangular then X is the Cholesky factor of A up to a column scaling bydiag(±1).The inverse of the Pascal matrix was found explicitly by Cohen [230, 1975]:the (i,j) element isThe Pascal matrix can be made singular simply by subtracting 1 from the(n,n) element.
To see this, note that26.4 T HE P ASCAL M ATRIX523This perturbation, ∆P = –enenT, is far from being the smallest one that makesP singular, which is ∆Popt =where λn is the smallest eigenvalue ofP and vn is the corresponding unit eigenvector, foris of order (n!)2 /(2 n)! ~as we saw above.A more subtle property of the Pascal matrix is that it is totally positive. Karlin [644, 1968, p. 137] shows that the matrix with elements( i, j = 0,1, .
. .) is totally positive; the Pascal matrix is a submatrix of this oneand hence is also totally positive. From the total positivity it follows that thePascal matrix has distinct eigenvalues, which (as we already know from thepositive definiteness) are real and positive, and that its ith eigenvector (corresponding to the ith largest eigenvalue) has exactly i – 1 sign changes [414,1959, Thin.
13, p. 105].T = pascal (n, 2) is obtained by rotating S clockwise through 90 degreesand multiplying by — 1 if n is even:>> T = pascal(6, 2)T=-15-1010-51-14-64-10-13-3100-12-1000-110000-100000It has the surprising property that it is a cube root of the identity, a propertynoted by Turnbull [1028, 1929, p. 332]:>> T*Tans =o0000-100001-1000-12-1001-33-10-14-64-11-510-lo5-1010001000000000>> T*T*Tans =100524A G ALLERYOFT EST M ATRICESFigure 26.1. spy(rem(pascal(32),2)).000000000100010001Finally, we note that it is trivial to plot an approximation to the Sierpinskigasket [824, 1992, §2.2], [478,1992] in MATLAB: simply type spy(rem(pascal (n),2)). See Figure 26.1. The picture produced by this command is incorrect forlarge n, however, because the elements of pascal(n) become too large to beexactly representable in floating point arithmetic.26.5. Tridiagonal Toeplitz MatricesA tridiagonal Toeplitz matrix has the formSuch matrices arise, for example, when discretizing partial differential equations or boundary value problems for ordinary differential equations.
The26.6 C OMPANION M ATRICES525eigenvalues are known explicitly [884, 1947], [885, 1952], [1006, 1977, pp. 155–156] :d + 2(ce)1/2 cos(kπ/(n + 1)),k = 1:n.The eigenvalues are also known for certain variations of the symmetric matrixT n ( c, d, c) in which the (1,1) and (n, n) elements are modified; see Gregoryand Karney [482, 1969].The matrix Tn (–1, 2, –1) is minus the well-known second difference matrix, which arises in applying central differences to a second derivative operator. Its inverse has (i, j) element –i(n – j + 1)/(n + 1) (cf.
Theorem 14.9).The condition number satisfies k2 (Tn ) ~ (4/π 2 )n 2 .One interesting property of Tn(c, d, e) is that the diagonals of its LU factorization converge as nwhen Tn is symmetric and diagonally dominant,and this allows some savings in the computation of the LU factorization, asshown by Malcolm and Palmer [725, 1974]. Similar properties hold for cyclicreduction; see Bondeli and Gander [134, 1994].26.6. Companion MatricesThe companion matrix associated with the characteristic polynomialof Ais the matrixThe Test Matrix Toolbox function compan computes C, via the call C =compan (A), It is easy to check that C has the same characteristic polynomial as A, and that if λ is an eigenvalue of C thenishas rank at least n – 1 for anya corresponding eigenvector. Since Cλ, C is nonderogatory, that is, in the Jordan form no eigenvalue appears inmore than one Jordan block.
It follows that A is similar to C only if A isnonderogatory.There are no explicit formulae for the eigenvalues of C, but, perhaps surprisingly, the singular values have simple representations, as found by Kenneyand Laub [651, 1988] (see also Kittaneh [660, 1995]):A G ALLERY526OFT EST M ATRICESFigure 26.2. Pseudospectra of compan(A).whereThe compan function is a useful means for generating new test matrices, compan(A) is a nonnormal matrix with the samefrom old.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
