# This code is stored in the file # # matrix.ssc # #------------------------------- # Add and subtract matrices #------------------------------- > a<-matrix(c(3, 6, 2, 1),2,2,byrow=T) > a [,1] [,2] [1,] 3 6 [2,] 2 1 > aa<-matrix(c(3, 6, 2, 1),2,2,byrow=F) > aa [,1] [,2] [1,] 3 2 [2,] 6 1 > b<-matrix(c(7, -4, -3, 2),2,2,byrow=T) > b [,1] [,2] [1,] 7 -4 [2,] -3 2 > a+b [,1] [,2] [1,] 10 2 [2,] -1 3 > a-b [,1] [,2] [1,] -4 10 [2,] 5 -1 #--------------------------------- # Multiplication by a scalar #--------------------------------- > c<-matrix(c(2, -1, 3, 0, 4, -2),2,3,byrow=T) > c [,1] [,2] [,3] [1,] 2 -1 3 [2,] 0 4 -2 > d<-2*c > d [,1] [,2] [,3] [1,] 4 -2 6 [2,] 0 8 -4 #------------------------------- # Transpose of a matrix #------------------------------- > ct <- t(c) > ct [,1] [,2] [1,] 2 0 [2,] -1 4 [3,] 3 -2 #------------------------------- # Matrix multiplication #------------------------------- > a<- matrix(c(3, 0, -2, 1, -1, 4), 2,3,byrow=T) > a [,1] [,2] [,3] [1,] 3 0 -2 [2,] 1 -1 4 > b<- matrix(c(1,1,1,2,1,3), 3,2,byrow=T) > b [,1] [,2] [1,] 1 1 [2,] 1 2 [3,] 1 3 > c<-a%*%b > c [,1] [,2] [1,] 1 -3 [2,] 4 11 #---------------------- # Inner product #---------------------- > x<-c(1,7,-6,4) > y<-c(2,-2,1,5) > x [1] 1 7 -6 4 > y [1] 2 -2 1 5 > t(x)%*%y [,1] [1,] 2 > x%*%y [,1] [1,] 2 > crossprod(x,y) [,1] [1,] 2 #---------------------- # Length of a vector #---------------------- > ynorm<-sqrt(crossprod(y,y)) > ynorm [,1] [1,] 5.830952 #---------------------------------- # Number of elements in a vector #---------------------------------- > length(y) [1] 4 #-------------------------------- # Elementwise multiplication #-------------------------------- > a<-matrix(c(1,2,3,4),2,2,byrow=T) > a [,1] [,2] [1,] 1 2 [2,] 3 4 > b<-matrix(c(3,-1,0,5),2,2,byrow=T) > b [,1] [,2] [1,] 3 -1 [2,] 0 5 > a*b [,1] [,2] [1,] 3 -2 [2,] 0 20 # ------------------------- # Kronecker Product) # ------------------------- > a <- matrix(c(2,4,0,-2,3,-1),ncol=2,byrow=T) > a [,1] [,2] [1,] 2 4 [2,] 0 -2 [3,] 3 -1 > b <- matrix(c(5,3,2,1),2,2,byrow=T) > b [,1] [,2] [1,] 5 3 [2,] 2 1 > kronecker(a,b) [,1] [,2] [,3] [,4] [1,] 10 6 20 12 [2,] 4 2 8 4 [3,] 0 0 -10 -6 [4,] 0 0 -4 -2 [5,] 15 9 -5 -3 [6,] 6 3 -2 -1 #--------------------------------------- # What happens when the dimensions # of the matrices or vectors are # not appropriate for the operation #--------------------------------------- > a<-matrix(c(1, 1, 1, 2), 2, 2, byrow=T) > b<-matrix(c(3, 0, -2, 1, -1, 4), 2, 3, byrow=T) > a [,1] [,2] [1,] 1 1 [2,] 1 2 > b [,1] [,2] [,3] [1,] 3 0 -2 [2,] 1 -1 4 > a+b Error in a + b: Dimension attributes do not match > b+a Error in b + a: Dimension attributes do not match > a%*%b [,1] [,2] [,3] [1,] 4 -1 2 [2,] 5 -2 6 > b%*%a Error in "%*%.default"(b, a): Number of columns of x should be the same as number of rows of y > a*b Error in a * b: Dimension attributes do not match # This sciprt is also stored in # # matrix.ssc #---------------------------- # Inverse of a matrix #---------------------------- > w<-matrix(c(1,2,3,4,5,6,7,8,10),3,3,byrow=T) > w [,1] [,2] [,3] [1,] 1 2 3 [2,] 4 5 6 [3,] 7 8 10 > winv<-solve(w) > winv [,1] [,2] [,3] [1,] -0.6666667 -1.333333 1 [2,] -0.6666667 3.666667 -2 [3,] 1.0000000 -2.000000 1 > w%*%winv [,1] [,2] [,3] [1,] 1.000000e+00 4.440892e-15 -2.664535e-15 [2,] 8.881784e-16 1.000000e+00 -8.881784e-16 [3,] 0.000000e+00 0.000000e+00 1.000000e+00 # ------------------------------ # Determinant of a matrix # ------------------------------ # Build your own function > determ<-function(M) Re(prod(eigen(M, only.values=T)$values)) > determ(w) [1] -3 # Another function (V&R, page 101) > absdet <- function(M) abs(prod(diag(qr(M)$qr))) > absdet(w) [1] 3 # Another example > x1 <- matrix(c(1,2,3,4,5,6,7,8,9),ncol=3,byrow=T) > x1 [,1] [,2] [,3] [1,] 1 2 3 [2,] 4 5 6 [3,] 7 8 9 > determ(x1) [1] 3.154999e-15 > absdet(x1) [1] 1.631688e-15 #-------------------------------------------- # Rank of a matrix: use the "qr" function # (V&R on p.101) #-------------------------------------------- > A <- matrix(c(1, 1, 1, + 2, 5, -1, + 0, 1, -1),3,3,byrow=T) > A [,1] [,2] [,3] [1,] 1 1 1 [2,] 2 5 -1 [3,] 0 1 -1 > qr(A)$rank [1] 2 # Another example > A <- matrix(c(1,1, 1, + 2,5,-1, + 0,1, 1),3,3,byrow=T) > A [,1] [,2] [,3] [1,] 1 1 1 [2,] 2 5 -1 [3,] 0 1 1 > qr(A)$rank [1] 3 # Another example > X <- matrix(c(1,1,0,0, + 1,1,0,0, + 1,0,1,0, + 1,0,1,0, + 1,0,0,1, + 1,0,0,1),ncol=4,byrow=T) > X [,1] [,2] [,3] [,4] [1,] 1 1 0 0 [2,] 1 1 0 0 [3,] 1 0 1 0 [4,] 1 0 1 0 [5,] 1 0 0 1 [6,] 1 0 0 1 > qr(X)$rank [1] 3 # Note that the sum of squares # and crossproducts matrix has # the same rank as X > XtX <- t(X)%*%X > XtX [,1] [,2] [,3] [,4] [1,] 6 2 2 2 [2,] 2 2 0 0 [3,] 2 0 2 0 [4,] 2 0 0 2 >> qr(XtX)$rank [1] 3 # This is a square symmetric matrix # but the inverse does not exist > solve(XtX) Problem in solve.qr(a): apparently singular matrix Use traceback() to see the call stack # Note that the function "rank" in Splus # is related to sorting. It computes the # ranks of the elements of a vector. # (V&R on p.53) > rank(c(1.2, 5.1, 3.5, 9.8)) [1] 1 3 2 4 > # --------------------------------- > # Create an nxn identity matrix > # --------------------------------- > diag(rep(1,4)) [,1] [,2] [,3] [,4] [1,] 1 0 0 0 [2,] 0 1 0 0 [3,] 0 0 1 0 [4,] 0 0 0 1 #--------------------------- # Trace of a matrix #--------------------------- > w<-matrix(c(1,2,3,4,5,6,7,8,10),3,3,byrow=T) > w [,1] [,2] [,3] [1,] 1 2 3 [2,] 4 5 6 [3,] 7 8 10 > tr<-sum(diag(w)) > tr [1] 16 #---------------------------------- # Compute row sums or column sums #---------------------------------- > sum(w) [1] 46 > apply(w,1,sum) [1] 6 15 25 > apply(w,2,sum) [1] 12 15 19 > apply(w,1,prod) [1] 6 120 560 > apply(w,1,mean) [1] 2.000000 5.000000 8.333333 > apply(w,1,var) [1] 1.000000 1.000000 2.333333 # ---------------------------------- # Eigenvalues & Eigenvectors # ---------------------------------- > A <- matrix(c(1.96,.72,.72,1.54),2,2,byrow=T) > A [,1] [,2] [1,] 1.96 0.72 [2,] 0.72 1.54 > EA <- eigen(A) > EA $values: [1] 2.5 1.0 $vectors: [,1] [,2] [1,] -0.8 0.6 [2,] -0.6 -0.8 > EA$values [1] 2.5 1.0 > eigen(A)$values [1] 2.5 1.0 # --------------------------------- # Singular Value Decomposition # --------------------------------- > A <- matrix(c(2,0,1,1,0,2,1,1,1,1,1,1), + ncol=4,byrow=T) > A [,1] [,2] [,3] [,4] [1,] 2 0 1 1 [2,] 0 2 1 1 [3,] 1 1 1 1 > svdA <- svd(A) > svdA $d: [1] 3.464102 2.000000 0.000000 $v: [,1] [,2] [,3] [1,] -0.5 7.071068e-001 0.5 [2,] -0.5 -7.071068e-001 0.5 [3,] -0.5 -1.226910e-016 -0.5 [4,] -0.5 -9.065285e-017 -0.5 $u: [,1] [,2] [,3] [1,] -0.5773503 7.071068e-001 -0.4082483 [2,] -0.5773503 -7.071068e-001 -0.4082483 [3,] -0.5773503 -7.597547e-017 0.8164966 > svdA$u %*% t(svdA$u) [,1] [,2] [,3] [1,] 1.000000e+000 9.310586e-018 -3.089976e-018 [2,] 9.310586e-018 1.000000e+000 3.244475e-017 [3,] -3.089976e-018 3.244475e-017 1.000000e+000 > t(svdA$v) %*% svdA$v [,1] [,2] [,3] [1,] 1.000000e+000 5.116079e-017 -2.775558e-017 [2,] 5.116079e-017 1.000000e+000 -3.405750e-017 [3,] -2.775558e-017 -3.405750e-017 1.000000e+000 > svdA$u%*%diag(svdA$d)%*%t(svdA$v) [,1] [,2] [,3] [,4] [1,] 2.000000e+000 -1.557456e-016 1 1 [2,] -9.074772e-017 2.000000e+000 1 1 [3,] 1.000000e+000 2.000000e+000 1 1 > diag(svdA$d) %*% diag(svdA$d) [,1] [,2] [,3] [1,] 12 0 0 [2,] 0 4 0 [3,] 0 0 0 > eigen(A%*%t(A))$values [1] 1.200000e+001 4.000000e+000 4.440892e-016 > eigen(t(A)%*%A)$values [1] 1.200000e+001 4.000000e+000 -2.167786e-016 -3.238078e-015 # An example where the singular values # are the eigenvalues > A <- matrix(c(1.96,.72,.72,1.54),2,2,byrow=T) > A [,1] [,2] [1,] 1.96 0.72 [2,] 0.72 1.54 > > svdA <- svd(A) > svdA $d: [1] 2.5 1.0 $v: [,1] [,2] [1,] -0.8 -0.6 [2,] -0.6 0.8 $u: [,1] [,2] [1,] -0.8 -0.6 [2,] -0.6 0.8 #---------------------------------------- # Trace and determinant of a matrix #---------------------------------------- > A <- matrix(c(1,1, 1, + 2,5,-1, + 0,1, 1),3,3,byrow=T) > A [,1] [,2] [,3] [1,] 1 1 1 [2,] 2 5 -1 [3,] 0 1 1 > traceA <- sum(diag(A)) > traceA [1] 7 > eigenA <- eigen(A) > eigenA $values: [1] 5.336912+0.0000000i 0.831544-0.6578603i 0.831544+0.6578603i $vectors: [,1] [,2] [,3] [1,] -0.2852888+0i 1.7077352+0.3055786i 1.7077352-0.3055786i [2,] -1.0054394+0i -0.7100770-0.2551929i -0.7100770+0.2551929i [3,] -0.2318330+0i 0.6234268-0.9197348i 0.6234268+0.9197348i > traceA <- sum(eigenA$values) > traceA [1] 7+0i > Re(traceA) [1] 7 > detA <- Re(prod(eigenA$values)) > detA [1] 6 # An example where the eigenvalues # are real numbers > A <- matrix(c(1,1, 1, + 2,5,-1, + 0, 1, -1),3,3,byrow=T) > A [,1] [,2] [,3] [1,] 1 1 1 [2,] 2 5 -1 [3,] 0 1 -1 > eigenA <- eigen(A) > eigenA $values: [1] 5.372281e+00 -3.722813e-01 -1.405092e-16 $vectors: [,1] [,2] [,3] [1,] -0.2608539 -3.269275 -0.8164966 [2,] -0.9858217 1.730136 0.4082483 [3,] -0.1547047 2.756228 0.4082483 > traceA <- sum(eigenA$values) > traceA [1] 5 > detA <- Re(prod(eigenA$values)) > detA [1] 2.810183e-16 #------------------------------------------- # Eigenvalues of a square symmetric matrix #------------------------------------------- > A<-matrix(c(4,2,-1,2,6,-4,-1,-4,9),3,3,byrow=T) > A > EA <- eigen(A) > EA $values: [1] 12.245772 4.433349 2.320879 $vectors: [,1] [,2] [,3] [1,] -0.2347350 -0.7321107 0.6394634 [2,] -0.5764345 -0.4248579 -0.6980108 [3,] 0.7827022 -0.5324563 -0.3222848 > SVDA <- svd(A) > SVDA $d: [1] 12.245772 4.433349 2.320879 $v: [,1] [,2] [,3] [1,] -0.2347350 0.7321107 -0.6394634 [2,] -0.5764345 0.4248579 0.6980108 [3,] 0.7827022 0.5324563 0.3222848 $u: [,1] [,2] [,3] [1,] -0.2347350 0.7321107 -0.6394634 [2,] -0.5764345 0.4248579 0.6980108 [3,] 0.7827022 0.5324563 0.3222848 #----------------------------------------- # An example of a square symmetric matrix # that is not positive definite #----------------------------------------- > W<-matrix(c(4,2,-1,2,6,-4,-1,-4,-9),3,3,byrow=T) > W [,1] [,2] [,3] [1,] 4 2 -1 [2,] 2 6 -4 [3,] -1 -4 -9 > EW <- eigen(W) > EW $values: [1] 8.151345 2.865783 -10.017128 $vectors: [,1] [,2] [,3] [1,] 0.4665008 -0.88381658 0.0352886 [2,] 0.8550024 0.46079428 0.2379907 [3,] -0.2266009 -0.08085101 0.9706262 > t(EW$vectors)%*%EW$vectors [,1] [,2] [,3] [1,] 1.000000e+000 -1.394453e-016 -1.517341e-016 [2,] -1.394453e-016 1.000000e+000 -5.759824e-017 [3,] -1.517341e-016 -5.759824e-017 1.000000e+000 > SVDW <- svd(W) > SVDW $d: [1] 10.017128 8.151345 2.865783 $v: [,1] [,2] [,3] [1,] -0.0352886 -0.4665008 -0.88381658 [2,] -0.2379907 -0.8550024 0.46079428 [3,] -0.9706262 0.2266009 -0.08085101 $u: [,1] [,2] [,3] [1,] 0.0352886 -0.4665008 -0.88381658 [2,] 0.2379907 -0.8550024 0.46079428 [3,] 0.9706262 0.2266009 -0.08085101 #-------------------------------------- # Inverse of a matrix #-------------------------------------- > A <- matrix(c(1.96,.72,.72,1.54),2,2,byrow=T) > Ainv <- solve(A) > Ainv [,1] [,2] [1,] 0.616 -0.288 [2,] -0.288 0.784 > A%*%Ainv [,1] [,2] [1,] 1.000000e+000 -1.638772e-016 [2,] 7.548758e-017 1.000000e+000 # Use the spectral decomposition # to compute the inverse of a matrix > Aev<-eigen(A)$vectors > Aeval<-eigen(A)$values > Ainv2<-Aev%*%solve(diag(Aeval))%*%t(Aev) > Ainv2 [,1] [,2] [1,] 0.616 -0.288 [2,] -0.288 0.784 #------------------------------------------ # Solutions to linear equations #------------------------------------------ > x<-c(1,1) > x [1] 1 1 > b<-solve(A,x) > b [1] 0.328 0.496 [2,] -0.288 0.784 > A%*%Ainv [,1] [,2] [1,] 1.000000e+000 -2.220446e-016 [2,] 5.551115e-017 1.000000e+000