install.packages("Matrix") library(Matrix) ## Matrix 1 A1 <- matrix(c(13,6,5, 6,12,6, 5,6,13), 3, 3) b1 <- matrix(c(35, 18, 43), ncol=1) ## Obtain Decompositions A1.inv <- solve(A1) ## B&M Section 6.2.3. pp. 169-170 A1.lu <- lu(A1) ## B&M Section 6.2.2. pp. 168-169 (Matrix package) A1.qr <- qr(A1) ## B&M Section 6.4.3. pp. 174-178 A1.chol <- chol(A1) ## B&M Section 6.4.2. pp. 173-174 A1.eigen <- eigen(A1) ## B&M Section 6.3. p. 171 A1.svd <- svd(A1) ## B&M Section 6.4.1. pp. 172-173 ## Extract and use Decompositions # "Basic inverse" and solution to Ax=b A1.inv A1 %*% A1.inv solve(A1, b1) A1.inv %*% b1 ## LU Decomposition and solution to Ax=b eluA1 <- expand(A1.lu) (L <- eluA1$L) (U <- eluA1$U) L %*% U (y <- solve(L,b1)) solve(U,y) ## QR Decomposition and solution to Ax=b (Q <- qr.Q(A1.qr)) (R <- qr.R(A1.qr)) Q %*% R solve(R, t(Q) %*% b1) qr.solve(R) %*% t(Q) ## Eigenvalue/Eigenvector Decomposition A1.eigen L <- diag(A1.eigen$values) Q <- A1.eigen$vectors Q %*% L %*% t(Q) ## Singular Value Decomposition A1.svd A1.svd$u %*% diag(A1.svd$d) %*% t(A1.svd$v) ## SVD - Direct calculation A1pA1.eigen <- eigen(t(A1) %*% A1) A1A1p.eigen <- eigen(A1 %*% t(A1)) U <- A1A1p.eigen$vectors V <- A1pA1.eigen$vectors D <- diag(sqrt(A1pA1.eigen$values)) U %*% D %*% t(V)