###################################################### # # R implementation of a PCA corresponding to exercises # 8.1, 8.2 and 8.3 of # http://www.dbs.ifi.lmu.de/cms/Maschinelles_Lernen_und_Data_Mining_12 # # programming language R available at # http://cran.r-project.org/ # # available in the CIP pool by # $ R # # This file was created by Marisa Petri (contact: # petri[aT}dbs.ifi.lmu.de) # ###################################################### # # computes the principal components of input matrix X # # X : a M x N matrix with N M-dimensional patterns # contained in the columns pca <- function(X) { # first: determine mean input vector mu <- apply(X, 1, mean) # define difference vectors (x_i - mu) Xd <- X - mu %*% t(rep(1,ncol(X))) # covariance matrix: 1/N Xd * Xd^t xcov <- 1/ncol(Xd) * (Xd %*% t(Xd)) # run eigenvalue decomposition evd <- eigen(xcov) ## Xd * Xd^t = U * D * U^t return (list(X=X, Xd=Xd, mu=mu, evd=evd, cov=xcov)) } p <- pca(t(num)) # # reduce the dimensionality of the input data matrix X to # the first r principal components # # r: the target number of principal components # Xd: the normalized input data matrix; xd_i = x_i - mu # U: eigenvectors of covariance matrix 1/N (Xd * Xd^t) # reduceDimensions <- function(r, Xd=p$Xd, U=p$evd$vectors) { V <- U[,1:r] # first r eigenvectors of U; M x r Xl <- t(V) %*% Xd # down-project X; r x N return(Xl) } # # reconstruction of input matrix X with r principal components # # r: the target number of principal components # X: the input data matrix (N column-wise M-dimensional entries) # Xd: the normalized input data matrix; xd_i = x_i - mu # mu: the mean vector of X (M dimensions) # U: eigenvectors of covariance matrix 1/N (Xd * Xd^t) # testReconstruct <- function(r, X=p$X, Xd=p$Xd, mu=p$mu, U=p$evd$vectors) { stopifnot(r > 0 && r <= nrow(X) && all(dim(X) == dim(Xd)) && length(mu) == nrow(X)) # etc. ... Xl <- reduceDimensions(r, Xd, U) # down-project X; r x N # now: reconstruction V <- U[,1:r] # first r eigenvectors of U; M x r # remember: Xd = X - mu * 1 # => reconstruct matrix X as # mu * 1 + V * V^t * Xd = mu * 1 + V * Xl rec <- mu %*% t(rep(1,ncol(Xd))) + V %*% Xl dffs <- c() # compute errors for (i in 1:ncol(X)) { dffs <- c(dffs, sum((X[,i] - rec[,i])^2)) } return(list(rec=rec, total_err=sum(abs(X - rec)), dffs=dffs)) } # exercise 8-1: # toy dataset with one dominating principal component S <- rbind( c(1:20, 4:15, 3:13, 8:12, 12:16), c(-(1:20) + 20, -(4:15) + 21, -(3:13) + 19, 8:12 + 6, 12:16 - 12) ) previousPar <- par(mfrow=c(1,2)) plot(S[1,], S[2,], main="original dataset in 2D", col=c(rep(1,43), rep(2,10)), lwd=2) p <- pca(S) rec <- testReconstruct(1)$rec plot(rec[1,], rec[2,], main="dataset reconstructed from \n1st principal component", col=c(rep(1,43), rep(2,10)), lwd=2) par(previousPar) # generate outlier So <- S; So[2,1] <- -50 previousPar <- par(mfrow=c(1,2)) plot(So[1,], So[2,], main="original dataset in 2D with outlier", col=c(rep(1,43), rep(2,10)), lwd=2) p <- pca(So) recO <- testReconstruct(1)$rec plot(recO[1,], recO[2,], main="dataset reconstructed from \n1st principal component", xlim=range(So[1,]), ylim=range(So[2,]), col=c(rep(1,43), rep(2,10)), cex=3) points(rec[1,], rec[2,], pch=20, col=c(rep(1,43), rep(2,10)), cex=3) legend("bottomleft", c("outlier reconstruction", "reconstruction without outliers"), pt.cex=3, pch=c(1,20)) par(previousPar) # exercise 8-2: X <- rbind(c(1:3,5:7), rep(c(0,6), each=3)) p <- pca(X) p # alternative way of eigenvalue computation # using det(Sigma - lambda I) = 0 ... (41/3 + sqrt((41/3)^2 - 24)) / 2 # lambda_1 (41/3 - sqrt((41/3)^2 - 24)) / 2 # lambda_2 # eigenvectors must then be solved from # X * u_j = lambda_j * u_j u <- p$evd$vectors d <- diag(p$evd$values) u %*% d %*% t(u) p$cov # dataset reduced to 1 principal component: reduceDimensions(1, p$Xd, p$evd$vectors) previousPar <- par(mfrow=c(1,2)) plot(X[1,], X[2,], col=1:6, lwd=2) text(p$mu[1], p$mu[2], expression(mu), cex=3) nU <- p$evd$vectors nU[,1] <- nU[,1] / sqrt(sum(nU[,1]^2)) nU[,2] <- nU[,2] / sqrt(sum(nU[,2]^2)) arrows(p$mu[1] - nU[1,1] * sqrt(p$evd$values[1]), p$mu[2] - nU[2,1] * sqrt(p$evd$values[1]), p$mu[1] + nU[1,1] * sqrt(p$evd$values[1]), p$mu[2] + nU[2,1] * sqrt(p$evd$values[1]), code=3) arrows(p$mu[1] - nU[1,2] * sqrt(p$evd$values[2]), p$mu[2] - nU[2,2] * sqrt(p$evd$values[2]), p$mu[1] + nU[1,2] * sqrt(p$evd$values[2]), p$mu[2] + nU[2,2] * sqrt(p$evd$values[2]), code=3) rec <- testReconstruct(1)$rec plot(rec[1,], rec[2,], col=1:6, lwd=2) par(previousPar) # exercise 8-3: num <- matrix(c( 0,1,0,0,1,0,0,1,0,0,1,0,0,1,0, # 1 1,1,1,0,0,1,1,1,1,1,0,0,1,1,1, # 2 1,1,1,0,0,1,0,1,1,0,0,1,1,1,1, # 3 1,0,1,1,0,1,1,1,1,0,0,1,0,0,1, # 4 1,1,1,1,0,0,1,1,1,0,0,1,1,1,1, # 5 1,1,1,1,0,0,1,1,1,1,0,1,1,1,1, # 6 1,1,1,0,0,1,0,0,1,0,0,1,0,0,1, # 7 1,1,1,1,0,1,1,1,1,1,0,1,1,1,1, # 8 1,1,1,1,0,1,1,1,1,0,0,1,1,1,1, # 9 1,1,1,1,0,1,1,0,1,1,0,1,1,1,1 # 0 ), ncol=15, byrow=TRUE) # visualize numbers: plotNumber <- function(number, x=num, dx=3, dy=5) { stopifnot(length(x[number,]) == dx * dy) # test dimensions # the matrix must be rotated by 90°; transpose causes the # image to be flipped, thus the flip transformation [dy:1,] image(t(matrix(x[number,],ncol=dx,byrow=T)[dy:1,])) } # visualize all in an array of images plotNumbers <- function(x=num, dx=3, dy=5) { previousPar <- par(mfrow=c(5,2), # configure the array mar = c(1, 1, 1, 1) + 0.1) # cut space around the edges for (i in 1:10) plotNumber(i, x, dx, dy) par(previousPar) # reset image device to default values } plotNumbers() # run PCA p <- pca(t(num)) # visualize eigenfaces and mean pixel values previousPar <- par(mfrow=c(4,4), # configure the array mar = c(1, 1, 3, 1) + 0.1) # cut space around the edges plotNumber(1, matrix(p$mu, nrow=1)) title("means") for (i in 1:ncol(num)) { plotNumber(i, t(p$evd$vectors)) title(paste("eigen value",i,":", round(p$evd$values[i],5))) } par(previousPar) # reset image device to default values # reconstruction of the input matrix: testReconstruct(15) # no (non-rounding-caused) error testReconstruct(8) # neither testReconstruct(7) # errors visisble for (i in ncol(num):1) { res <- testReconstruct(i) cat(sprintf("err d = %2d: %2.4f, avg(DFFS) = %2.4f \n",i,res$total_err, mean(res$dffs))) } # why? eigenvales: p$evd$values # visual analysis: visualizeReconstruct <- function(r) { tr <- testReconstruct(r) plotNumbers(t(tr$rec)) cat(sprintf("total absolute reconstruction error = %.2f\ndffs = %s\navg dffs = %.2f", tr$err, paste(c(1:9,0), ":", tr$dffs, collapse=","), mean(tr$dffs))) } visualizeReconstruct(15) visualizeReconstruct(8) visualizeReconstruct(7) visualizeReconstruct(5) visualizeReconstruct(1)