{VERSION 3 0 "IBM INTEL NT" "3.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 }{CSTYLE "" -1 23 "Courier" 1 10 0 0 0 0 0 0 0 0 0 0 3 0 0 } {CSTYLE "Help Normal" -1 30 "Times" 1 12 0 0 0 0 0 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 256 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 257 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 258 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 259 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 } {CSTYLE "" -1 260 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 261 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 262 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 263 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 264 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 265 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 266 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 267 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 } {CSTYLE "" -1 268 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 269 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Heading 1" 0 3 1 {CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 1 0 0 0 0 0 0 0 }1 0 0 0 8 4 0 0 0 0 0 0 -1 0 }{PSTYLE "Heading 2" 3 4 1 {CSTYLE "" -1 -1 "" 1 14 0 0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 8 2 0 0 0 0 0 0 -1 0 }{PSTYLE "Bullet Item" 0 15 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 3 3 0 0 0 0 0 0 15 2 }} {SECT 0 {EXCHG {PARA 0 "" 0 "" {TEXT -1 89 "MATH 2010 \+ Handout 3 - Condition Number of a Matrix" }} {PARA 0 "" 0 "" {TEXT -1 9 "Fall 2001" }}{PARA 0 "" 0 "" {TEXT -1 17 " Dr. R. Ablamowicz" }}{PARA 0 "" 0 "" {TEXT -1 98 "\nIn Section 2.3 I a ssigned problems 37, 38, 40, and 41, pp. 124-125. These problems deal \+ with the " }{TEXT 256 16 "condition number" }{TEXT -1 310 " of a matri x and the issue of sensitivity of the solution of a system of equation s to changes in the number of places of accuracy in the input. Before \+ we define the condition number, we solve problem 37.\n\nReference:\n\n [1] \"Linear Algebra and Its Applications\", 2nd edition, David C. La y, Addison-Wesley, 2000.\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "restart:with(linalg):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 258 30 "Pro blem 37 page 124 (see [1])." }{TEXT -1 10 "\n\nPart (a)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "eq1:=4.5*x1+3.1*x2=19.249;\neq2:=1. 6*x1+1.1*x2=6.843;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 64 "The simples t way to solve this system is to use command 'solve':" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "sol1:=solve(\{eq1,eq2\},\{x1,x2\}); " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 179 "Let's remember the coefficie nt matrix of this system as matrix A1. We will use it below after we d iscuss the condition number. Observe how small the determinant of this matrix is!" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "A1:=genmatri x(\{eq1,eq2\},[x1,x2]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "d et(A1);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "rank(A1);" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 5 "Now, " }{TEXT 269 104 "we change eq uations using only two decimal places of accuracy in the right-hand-si de of the above system" }{TEXT -1 1 ":" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "eq11:=4.5*x1+3.1*x2=19.25;\neq22:=1.6*x1+1.1*x2=6.84; " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 35 "and we solve this new system \+ again:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 86 "A2:=genmatrix(\{e q1,eq2\},[x1,x2]);\ndet(A2);\nrank(A2);\nsol2:=solve(\{eq11,eq22\},\{x 1,x2\});" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 199 "Observe, how differe nt now solution sol1 is from the earlier solution displayed in sol2. A nd this is only because we have used two instead of three decimal plac es of accuracy in the right hand side!\n" }}{PARA 0 "" 0 "" {TEXT -1 287 "Thus, the solution to this last system is a very poor approximati on to the exact solution of the first system even though the entries i n the right-hand-side of the last system differ less than 0.05% from t he entries in the right-hand-side of the original system. Where is the problem? \n\n" }{TEXT 257 20 "The condition number" }{TEXT -1 1 "\n" }}{PARA 0 "" 0 "" {TEXT -1 48 "The following has been taken from [JH] \+ page 114." }}{PARA 0 "" 0 "" {TEXT -1 5 "\nThe " }{TEXT 260 17 "condit ion number " }{TEXT -1 215 "k of an invertible matrix A is defined to \+ be the product of \n\nk = norm(A)*norm(A^(-1)), \n\nwhere the matrix n orm(A) is defined by max(norm(AX)/norm(X)) and the maximum is taken ov er all nonzero column vectors X. The " }{TEXT 261 11 "vector norm" } {TEXT -1 157 " norm(X) is usually the Euclidean length, but in LINPACK (a standard linear algebra package for numerical computations program med in FORTRAN) it is equal to " }{TEXT 262 52 "\"the sum of the absol ute values of all entries in X\"" }{TEXT -1 74 " since this is faster \+ to compute and since the matrix norm is then simply " }{TEXT 263 54 " \"the maximum of the vector norms of all columns of A\"." }{TEXT -1 460 "\n\nThe computation of the matrix norm of the inverse A^(-1) of A is not so simple. In fact we, don't even want to assume that A^(-1) e xists, since we want a measure that will help us distinguish singular \+ matrices, i.e., non-invertible, from nonsingular, i.e., invertible, on es. Note, however that \n\n1/norm(A^(-1)) = 1/max(norm(A^(-1)X)/norm(X )) = min(norm(AY)/norm(Y)) \n\nand that this last expression is define d even when A is singular. We therefore define the " }{TEXT 264 11 "r- condition" }{TEXT -1 65 " of any square matrix A to be \n\n(1/norm(A)) *min(norm(AY)/norm(Y))" }}{PARA 0 "" 0 "" {TEXT -1 301 "\nand we note \+ that it equals 1/k when A is nonsingular. The above minimum can be app roximated by choosing X so that norm(X)/norm(Y) is as small as possibl e, where Y is a solution of AY = X. For more information please consul t the references below such as [DMBS] or Maple's help page [Maple] sho wn next." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 185 "In short, the larger the condition number of a matrix, the closer this matrix is to being singular (non-invertible). Substantial errors may be produced when a condition number is large." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 12 "linalg[cond]" } {TEXT 30 31 " - condition number of a matrix" }}{PARA 4 "" 0 "usage" {TEXT -1 17 "Calling Sequence:" }}{PARA 0 "" 0 "" {TEXT -1 12 " co nd(A)" }}{PARA 0 "" 0 "" {TEXT -1 22 " cond(A, normname)" }}{PARA 4 "" 0 "" {TEXT -1 11 "Parameters:" }}{PARA 0 "" 0 "" {TEXT -1 5 " \+ " }{TEXT 23 11 "A - " }{TEXT -1 16 "a square matrix " }}{PARA 0 "" 0 "" {TEXT -1 5 " " }{TEXT 23 11 "normname - " }{TEXT -1 73 " (optional) matrix norm, must be one of: 1, 2, 'infinity', or 'frobeniu s'." }}}{SECT 0 {PARA 4 "" 0 "info" {TEXT -1 12 "Description:" }} {PARA 15 "" 0 "" {TEXT -1 222 "The function cond computes the ``standa rd'' matrix condition number, defined as norm(A) * norm(inverse(A)). T he matrix norm is the default employed by the linalg[norm] function, n amely the infinity norm (maximum row sum). " }}{PARA 15 "" 0 "" {TEXT -1 176 "More generally, cond(A, normname) computes norm(A, normname) * norm(inverse(A), normname). This is the same measure, but using the s pecified norm instead of the infinity norm. " }}{PARA 15 "" 0 "" {TEXT -1 85 "The command with(linalg,cond) allows the use of the abbre viated form of this command." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "A:=matrix(3,3,[1,2,-3,4,5,6,7,8,-9]);" }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 100 "'1'-norm of a matrix is the maximum of the column sum \+ of magnitudes, that is used in LINPACK, hence:" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 10 "norm(A,1);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 108 "The '2'-norm of a matrix is the square root of the maximum eig envalue of the matrix A * transpose(A), hence:" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 17 "evalf(norm(A,2));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 137 "For a positive integer k, the k-norm of a vector is the \+ k-th root of the sum of the magnitudes of each element raised to the k -th power. " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "X:=vector([1 ,2,-3]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 34 "hence the 1-norm of a vector, , is" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "norm(X,1); #used in LINPACK" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 31 "while the 2- norm of a vector is" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "norm (X,2);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 241 "For vectors, the infin ity norm is the maximum magnitude (absolute values) of all elements. T he infinity norm of a matrix is the maximum row sum, where the row sum is the sum of the magnitudes (absolute values) of the elements in a g iven row:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "norm(X,infinit y);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "norm(A,infinity);" } }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 140 "The frobenius norm of a matrix \+ or vector is defined to be the square root of the sum of the squares o f the magnitudes of each element. Thus," }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "norm(X,frobenius);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "norm(A,frobenius);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 224 "Thus, by default, Maple computes the condition number using th e infinity norm. In order to make it use the 1-norm used in LINPACK, w e need to tell it to use the 1-norm. For example, the condition number of matrix A above is:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "c ond(A,1);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 82 "while the condition \+ number of an identity matrix is 1 in norms 1, 2, and infinity:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "B:=diag(1,1,1);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "cond(B,1),cond(B,2),cond(B,infinity ),cond(B,frobenius);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 74 "Let's tak e a look at the coefficient matrix A1 from problem 37 done above:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "cond(A1,1);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 65 "Notice that the determinant of A1 is rath er small, close to zero:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 " det(A1);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 75 "On the other hand, Ma ple still produces a correct value for the rank of A1:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "rank(A1);" }}}{EXCHG {PARA 0 "" 0 " " {TEXT 267 33 "Problem 38 on page 125 (see [1]):" }}}{EXCHG {PARA 0 " > " 0 "" {MPLTEXT 1 0 105 "M:=matrix(4,4,[7,-6,-4,1,\n - 5,1,0,-2,\n 10,11,7,-3,\n 19,9,7,1]);" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "det(M);" }}}{EXCHG {PARA 0 " " 0 "" {TEXT -1 48 "Notice that the condition number for M is large:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "cond(M,1);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "x:=matrix(4,1,[1.56,0.67,-10.34,5.7 8]); #random vector x" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 59 "b: =evalm(M &* x);b:=convert(b,vector);x:=convert(x,vector);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 88 "Now, we are supposed to solve system Ax = b for x, where A=M, and call the solution x1: " }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 12 "A:=evalm(M);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "sys:=geneqns(A,[x1,x2,x3,x4],b);" }}}{EXCHG {PARA 0 " > " 0 "" {MPLTEXT 1 0 30 "sol:=solve(sys,\{x1,x2,x3,x4\});" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "x1:=subs(sol,vector([x1,x2,x3,x4])) ;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "x=evalm(x);" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 151 "Thus, Maple has found exactly the same vector x, that is, accuracy is 100%. Hoever, it is likely that a nother software would not produce exact results." }}}{EXCHG {PARA 0 " " 0 "" {TEXT 266 30 "Problem 40 page 125 (see [1])." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "A:=hilbert(5);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "rank(A);det(A);evalf(det(A));cond(A);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "B:=diag(1,1,1,1,1);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "AB:=augment(A,B);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "C:=gaussjord(AB);" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 30 "Ainv:=submatrix(C,1..5,6..10);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "evalm(A &* Ainv);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 144 "Thus, Maple can handle Hilbert matrices exactl y. Let's repeat this problem but take, instead of fractional entries, \+ floating point entries in A:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "AB:=augment(map(evalf,A),B);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "C:=gaussjord(AB);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "Ainv:=submatrix(C,1..5,6..10);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 109 "When we compare the above Ainv with the one previou sly obtained we find that now accuracy is just two digits!" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 268 20 "Problem 40 page 125:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "A:=hilbert(12);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "ev alf(cond(A));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "Ainv:=inve rse(A):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "evalm(Ainv &* A) ;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "Ainv2:=inverse(map(eva lf,A)):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "norm(Ainv-Ainv2, 1);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "evalm(Ainv2 &* A);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 " " {TEXT 259 10 "References" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 194 "[JH] Jepsen, Charles H., and Herman, Eugene A.: \+ \"MAX: MAtriX Algebra Calculator - Linear Algebra Problems for Compute r Solution\", Brooks/Cole Publishing Company, Pacific Grove, Californi a, 1988." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 113 "[Maple] Maple V Rel. 5.1 Help page on condition number of a matri x [accessed via ?condition], Waterloo Maple Inc." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 65 "[DMBS] Dongarra, J.J., C. B. Moler, J.R. Bunch, and G.W. Stewart: " }{TEXT 265 20 "LINPACK Users ' Guide" }{TEXT -1 13 ", SIAM, 1979." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}}{MARK "0 0 0" 61 }{VIEWOPTS 1 1 0 1 1 1803 }