reconstruct.svd<-function(s,k,bound=F) { if (k==1) o=outer(s$u[,1]*s$d[1],s$v[,1]) else o=s$u[,1:k]%*%diag(s$d[1:k])%*%t(s$v[,1:k]) if (bound) { o[o<0]=0 o[o>1]=1 } o } a=matrix(1:20,5,4) s=svd(a) s$u[,1:2]%*%diag(s$d[1:2])%*%t(s$v[,1:2]) reconstruct.svd(s,1) reconstruct.svd(s,2) #Chess board (a=matrix(floor((1.1249*(1:64))%%2),8,8)) s=svd(a) round(s$d,8) round(s$u,8) round(s$v,8) reconstruct.svd(s,1) reconstruct.svd(s,2) #PCA reconstruct.pca <- function(p,x,k) { if (k==1) o=outer(predict(p,x)[,1],p$rotation[,1:k]) else o=predict(p,x)[,1:k]%*%t(p$rotation[,1:k]) o=t(o) if (p$scale) o=o*p$scale o=o+p$center t(o) } x=matrix(rnorm(20000),10000,2)%*%matrix(rnorm(4),2,2) x[,1]=x[,1]+100 x[,2]=x[,2]+200 plot(x,asp=1) p=prcomp(x,scale=F) p$sdev p$rotation p$center plot(p$x,asp=1) plot(reconstruct.pca(p,x,1),asp=1) plot(reconstruct.pca(p,x,2),asp=1) #By svd x.centered=t(t(x)-apply(x,2,mean)) plot(x.centered) s=svd(x.centered) pr.comp=s$d/sqrt(dim(x)[1]) plot(s$u%*%diag(s$d),asp=1)