Here is a MAPLE BUG for the benefit of the MUG. It’s still there in MAPLE V Release 4, although I ran this under UNIX maple. I read the output into the program (or the other way around).
Notice the MAPLE 2-norm is actually the _SMALLEST_ singular value, not the largest. The
2-condition is the reciprocal of what it should be, which is the ratio of the largest over the
smallest singular value.
# This MAPLE program illustrates a bug in MAPLE that occurs ## when computing 2-norms of matrices. restart; with(linalg): Warning: new definition for norm Warning: new definition for trace A:=randmatrix(5,5); [ -85 -55 -37 -35 97 ] [ ] [ 50 79 56 49 63 ] [ ] A := [ 57 -59 45 -8 -93 ] [ ] [ 92 43 -62 77 66 ] [ ] [ 54 -5 99 -61 -50 ] sgm:=evalf(Svd(A)); sgm := [ 215.3159815, 192.8527895, 110.7852504, 66.73636029, 33.05625793 ] norm2A_1:=sgm[1]; norm2A_1 := 215.3159815 norm2A_2:=evalf(norm(A,2)); norm2A_2 := 33.05625794 cond2_1:=sgm[1]/sgm[5]; cond2_1 := 6.513622381 cond2_2:=evalf(cond(A,2)); cond2_2 := .1535244049
It is corrected with Maple 6. (U. Klein)
One of my students drew my attention to your remark on Maple’s behaviour with respect to the 2-norm of matrices.
Let me refer to the nice book GOLUB/ORTEGA, Scientific Computing - An Introduction with Parallel Computing, Academic Press, 1993. There I find a short and very well usable chapter on basic facts that are of interest in the field of numerical linear algebra. The subchapter 2.3 is dealing with norms.
The authors define the EUCLIDean matrix norm - similar to the way Maple does. See then subchapter 6.2 for norm based condition numbers. Again, GOLUB and ORTEGA give the definition Maple is also using.
Trying ?linalg,norm , you will get the information:
"The '2'-norm of a matrix is the square root of the maximum eigenvalue of the matrix A * transpose(A) . ".
You may use your example to reproduce Maple’s definition:
> restart:with(linalg): Warning, new definition for norm Warning, new definition for trace > A:=randmatrix(5,5); [-85 -55 -37 -35 97] [ ] [ 50 79 56 49 63] [ ] A := [ 57 -59 45 -8 -93] [ ] [ 92 43 -62 77 66] [ ] [ 54 -5 99 -61 -50] > Automatic:=norm(A,2); Automatic := sqrt(max(| RootOf(-102991333311217647241 2 + 130759387784582485 _Z - 36880268261270 _Z 3 4 5 + 3286110146 _Z - 101373 _Z + _Z ) |)) > ByHand:=sqrt(max(eigenvals(evalm(A &* transpose(A))))); ByHand := sqrt(max(RootOf(-102991333311217647241 2 + 130759387784582485 _Z - 36880268261270 _Z 3 4 5 + 3286110146 _Z - 101373 _Z + _Z ))) > evalf(Automatic)=evalf(ByHand); 33.05625794 = 33.05625794
Trying ?linalg,condition , you get the text:
"cond(A, normname) computes norm(A, normname) * norm(inverse(A), normname)."
Again, you may use your example to test Maple:
> evalf(cond(A,2))=evalf(norm(A,2)*norm(A^(-1),2)); .1535244049 = .1535244049
I use the same book, subchapter 2.2, to recall that B:=A*transpose(A) and C:=transpose(A)*A
both are positive definite (which means that they cannot have negative eigenvalues). Moreover,
both of them are symmetric, which means that det(B-lambda*I)=det(C-lambda*I) holds
and, thus, B and C do have the same eigenvalues (not: eigenvectors!). Hence, the
GOLUB/ORTEGA definition (using the spectral radius of C to determine the EUCLID
norm of A) and Maple’s definition (using the spectral radius of B to do the same) coincide
and reflect what can be found elsewhere as well.
> ByHandTranspose:=sqrt(max(eigenvals(evalm(transpose(A)&*A)))); ByHandTranspose := sqrt(max(RootOf(-102991333311217647241 2 + 130759387784582485 _Z - 36880268261270 _Z 3 4 5 + 3286110146 _Z - 101373 _Z + _Z ))) > evalf(Automatic),evalf(ByHand),evalf(ByHandTranspose); 33.05625794, 33.05625794, 33.05625794
Taking all this into account, I am sure that you will agree, that Maple is quite good in calculating the EUCLIDean norm or the EUCLIDean condition number of a matrix.
I think there is no bug.
Nevertheless, it is also very common to use condition numbers different to that seen above. The definition you mention reminds me of the TODD number (R. TODD, The Condition of Certain Matrices, Quart. J. Mech. and Appl. Math. 2 (1949) 469-472; Math. Rev. 11, 619) which is given by the quotient of the largest of the absolute vaues of the eigenvalues of A over the smallest of the absolute vaues of the eigenvalues of A. This is, in other words, the spectral radius of A, multiplied by the spectral radius of the inverse of A.
In this case, your example yields:
> num:=max(op(map(abs,[allvalues(eigenvalues(A))]))); num := 125.5772260 > den:=min(op(map(abs,[allvalues(eigenvalues(A))]))); den := 78.82169898 > T:=num/den; T := 1.593180909 > fac:=max(op(map(abs,[allvalues(eigenvalues(A^(-1)))]))); fac := .01268686178 > T=num*fac; 1.593180909 = 1.593180909