Guide to do SVD and PCA in R.
First some good documents about PCA and SVD:
http://www.math.ucsd.edu/~gptesler/283/slides/pca_f13-handout.pdf … about genetic application, VERY GOOD
http://www.cs.cmu.edu/~jimeng/papers/ICMLtutorial.pdf … It is a generic document about tensors on machine learning but has an initial chapeter very interesting about SVD and PCA
Haz clic para acceder a PCA-Tutorial-Intuition_jp.pdf
Imagine you have an matrix Y (mxn) with m rows and n columns, suppose m>=n. For example it can be a matrix of images with the rows having the vectorised pixel intesities of images or a matrix of genes with rows being the genes expressions for different samples. The columns have samples, so , for example, different columns correspond to different tissue samples or different images….
The SVD of Y is a matrix fatorization: Y= UDV’ , where the prime indicates transpose
In R SVD is obtained with:
s = svd(Y)
s$u gives the U
s$d gives the digonal vector of D
s$v … here we have the V, not the transposed V
Before applying svd the matrix Y have tobe detrended, that is, we have to substract the mean (through the sample dimension…)
so…previously…
Y= Y – rowMeans(Y)
Once this is done….the Principal Components are the Vs obtained in s$v after applying svd to the detrended Y
These can be obtained also using the function prcomp in R
p = prcomp(y)
and the principal components are in p$rotation…. try this plot to see that the columns of s$v and p$rotation are equal:
plot(s$v[,1],p$rotation[,1]) ,same for the other columns….
Before applying prcomp it is not necessary to do detrending.
THE PCs can be the s$v or the s$u depending on the problem and if we are working with Y or Y’. In general there is a confusion between what is the PCs and what is the eigenvector you take to project in the new low-dimension space, and there is confusion between documents about this.
In principle, the projection is done with the values obtained from svd (s$s or s$u) or their transposed. These values are usually called principal components, but principal components are called, in other papers, the above mentioned values multiplied by s$d, in these cases can be an additional factor of 1/m-1 to account for the factor to obtain the covariance matrix. See the first reference where it is presented quite clearly the differences.
IN GENERAL do YY’ and Y’Y and take the one with the lower dimension (this is usually the covariance matrix that makes sense), after doing that if you have done YY’ the PCs are the Us and if Y’Y then the PCs are the Vs
Here is some code in R you can use to see the differences between prcomp and svd:
u=s$u
d=s$d
v=s$v
View(p$rotation)
View(s$v)
View(abs(s$v-p$rotation))
sqrt(colSums(p$rotation^2)) # unitary vectors
View(Y %*% p$rotation -p$x) # p$x has the projections of Y on the PCs
Machine Learning, Maths and Physics 