Why does Andrew Ng prefer to use SVD and not EIG of covariance matrix to do PCA? I am studying PCA from Andrew Ng's Coursera course and other materials. In the Stanford NLP course cs224n's first assignment, and in the lecture video from Andrew Ng, they do singular value decomposition instead of eigenvector decomposition of covariance matrix, and Ng even says that SVD is numerically more stable than eigendecomposition.
From my understanding, for PCA we should do SVD of the data matrix of (m,n) size, not of the covariance matrix of (n,n) size. And eigenvector decomposition of covariance matrix.
Why do they do SVD of covariance matrix, not data matrix?
 A: For Python users, I'd like to point out that for symmetric matrices (like the covariance matrix), it is better to use numpy.linalg.eigh function instead of a general numpy.linalg.eig function.
eigh is 9-10 times faster than eig on my computer (regardless of matrix size) and has better accuracy (based on @Aksakal's accuracy test).
I am not convinced with the demonstration of the accuracy benefit of SVD with small eigenvalues. @Aksakal's test is 1-2 orders of magnitude more sensitive to random state than to the algorithm (try plotting all errors instead of reducing them to one absolute maximum). It means that small errors in the covariance matrix will have a greater effect on accuracy than the choice of an eigendecomposition algorithm. Also, this is not related to the main question, which is about PCA. The smallest components are ignored in PCA.
A similar argument can be made about numerical stability. If I have to use the covariance matrix method for PCA, I would decompose it with eigh instead of svd. If it fails (which has not been demonstrated here yet), then it is probably worth rethinking the problem that you are trying to solve before starting to look for a better algorithm.
A: Some great answers already have been given to your questions, so I won't add a lot of new stuff. But I tried (i) to base my answer on the knowledge you seem to have and (ii) to be as concise as possible. So you - or others in a similar situation - may find this answer helpful.
(Simple) Mathematical Explanation
SVD and the eigendecomposition are closely related. Let $X \in \mathbb{R}^{n \times p}$ be a real data matrix, so you may define its covariance matrix $C \in \mathbb{R}^{p\times p}$ as
\begin{equation}
C = \frac{1}{n} X^T X.
\end{equation}
1 | SVD of X
As you correctly stated, applying SVD on $X$ decomposes your original data in
\begin{equation*}
X = U S V^T
\end{equation*}
with $U \in \mathbb{R}^{n \times n}$, $V \in \mathbb{R}^{p \times p}$ being unitary, containing the (orthonormal) principal components and eigenvectors, respectively. The diagonal matrix $S \in \mathbb{R}^{n \times p}$ holds the singular values $s$. 
2 | Eigendecomposition of C
Since $C$ is hermitian its eigendecomposition yields eigenvectors given by the unitary matrix $V$ with corresponding real eigenvalues $\lambda$ as entries of a diagonal matrix $\Lambda \in \mathbb{R}^{p \times p}$:
\begin{equation}
CV = V\Lambda
\end{equation}
In this case we may calculate the principal components by projecting the eigenvectors on the original data PCs $= X V^T$. Note that these PCs are scaled by their corresponding eigenvalues and are thus correlated.
3 | SVD of C
In order to answer your questions, recall that we can factorize $C$ - since it is symmetric - via
\begin{equation}
C = V\Lambda V^T
\end{equation}
using its eigenvectors and eigenvalues.
Note, that this true also by just rearranging the equation from section (2). We can therefore calculate the eigenvectors of $X$ by applying the SVD to $C$. 
With just a bit more of effort we can now establish the relation between the singular values and the eigenvalues. Using the definition of the covariance we may as well write:
\begin{align}
C &= \frac{1}{n} ( U S V^T )^T ( U S V^T ) \\
&= \frac{1}{n} ( V S U^T U S V^T ) \\
&= \frac{1}{n} ( V S^2 V^T )
\end{align}
The last equation holds since $U$ is unitary, that is $U^T U = \mathbb{1}$. Now by simply comparing this result with that from above we find:
\begin{equation}
\frac{1}{n} ( V S^2 V^T ) = V \Lambda V^T \quad \Rightarrow \quad \lambda = \frac{s^2}{n} 
\end{equation}

Python, NumPy and Algorithms
Just a basic example to explore the different behaviour of numpy's


*

*linalg.svd(X)

*linalg.svd(C)

*linalg.eig(C)

*linalg.eigh(C)
shows that linalg.eig() reveals some (at least for me) unexpected behaviour. Calculating the matrix $V^TV=\mathbb{1}$ for all four cases, we can get a visual idea of the respective precision. It seems from the figure below, that linalg.eig() only provides a stable solution up to dimension $d = \text{rank}(C) = \text{min}(n,p)$.


# Create random data
n,p = [100,300]
X = np.random.randn(n,p)

# Covariane matrix
C = X.T @ X /n


# Create figure environment
fig = plt.figure(figsize=(14,5))
ax1 = fig.add_subplot(141)
ax2 = fig.add_subplot(142)
ax3 = fig.add_subplot(143)
ax4 = fig.add_subplot(144)

# 1. SVD on X
# ---------------------
U,s,VT = np.linalg.svd(X)
V = VT.T
ax1.imshow(V.T@V,cmap='Reds',vmin=0,vmax=1)

# 2. SVD on C
# ---------------------
V,eigenvalues,VT = np.linalg.svd(C)
ax2.imshow(V.T@V,cmap='Reds',vmin=0,vmax=1)
# 3. Eigendecomposition on C
# -> linalg.eig()
# ---------------------
eigenvalues,V = np.linalg.eig(C)
sortIdx = np.argsort(eigenvalues)[::-1]
V = V[:,sortIdx]
ax3.imshow((V.T@V).real,cmap='Reds',vmin=0,vmax=1)

# 4. Eigendecomposition on C
# -> linalg.eigh()
# ---------------------
eigenvalues,V = np.linalg.eigh(C)
sortIdx = np.argsort(eigenvalues)[::-1]
V = V[:,sortIdx]
ax4.imshow((V.T@V).real,cmap='Reds',vmin=0,vmax=1)

for a in [ax1,ax2,ax3,ax4]:
    a.set_xticks([])
    a.set_yticks([])
ax1.set_title('svd(X)')
ax2.set_title('svd(C)')
ax3.set_title('eig(C)')
ax4.set_title('eigh(C)')
fig.subplots_adjust(wspace=0,top=0.95)  
fig.suptitle('Eigendecomposition in NumPy')
plt.show()

A: amoeba already gave a good answer in the comments, but if you want a formal argument, here it goes.
The singular value decomposition of a matrix $A$ is $A=U\Sigma V^T$, where the columns of $V$ are eigenvectors of $A^TA$ and the diagonal entries of $\Sigma$ are the square roots of its eigenvalues, i.e. $\sigma_{ii}=\sqrt{\lambda_i(A^TA)}$.
As you know, the principal components are the orthogonal projections of your variables onto the space of the eigenvectors of the empirical covariance matrix $\frac{1}{n-1}A^TA$. The variance of the components is given by its eigenvalues, $\lambda_i(\frac{1}{n-1}A^TA)$.
Consider any square matrix $B$, $\alpha \in \mathbb R$ and a vector $v$ such that $Bv=\lambda v$. Then


*

*$B^kv=\lambda^kv$

*$\lambda(\alpha B) = \alpha\lambda( B)$


Let us define $S=\frac{1}{n-1}A^TA$.  The SVD of $S$ will compute the eigendecomposition of $S^TS=\frac{1}{(n-1)^2}A^TAA^TA$ to yield 


*

*the eigenvectors of $(A^TA)^TA^TA=A^TAA^TA$, which by property 1 are those of $A^TA$ 

*the square roots of the eigenvalues of $\frac{1}{(n-1)^2}A^TAA^TA$, which by property 2, then 1, then 2 again, are $\sqrt{\frac{1}{(n-1)^2} \lambda_i(A^TAA^TA)} = \sqrt{\frac{1}{(n-1)^2} \lambda_i^2(A^TA)} = \frac{1}{n-1}\lambda_i(A^TA) = \lambda_i(\frac{1}{n-1}A^TA)$.


Voilà!
Regarding the numerical stability, one would need to figure out what the employed alogrithms are. If you're up to it, I believe these are the LAPACK routines used by numpy:


*

*numpy.linalg.eig

*numpy.linalg.svd
Update: On the stability, the SVD implementation seems to be using a divide-and-conquer approach, while the eigendecomposition uses a plain QR algorithm. I cannot access some relevant SIAM papers from my institution (blame research cutbacks) but I found something that might support the assessment that the SVD routine is more stable. 
In 

Nakatsukasa, Yuji, and Nicholas J. Higham. "Stable and efficient spectral divide and conquer algorithms for the symmetric eigenvalue decomposition and the SVD." SIAM Journal on Scientific Computing 35.3 (2013): A1325-A1349.

they compare the stability of various eigenvalue algorithms, and it seems that the divide-and-conquer approach (they use the same one as numpy in one of the experiments!) is more stable than the QR algorithm. This, together with claims elsewhere that D&C methods are indeed more stable, supports Ng's choice.
A: To answer the last part of your question, "Why do they do SVD of covariance matrix, not data matrix?" I believe it is for performance and storage reasons. Typically, $m$ will be a very large number and even if $n$ is large, we would expect $m \gg n$.
Calculating the covariance matrix and then performing SVD on that is vastly quicker than calculating SVD on the full data matrix under these conditions, for the same result.
Even for fairly small values the performance gains are 
factors of thousands (milliseconds vs seconds). I ran a few tests on my machine to compare using Matlab:
 
That's just CPU time, but storage needs are just as, if not more, important. If you attempt SVD on a million by a thousand matrix in Matlab it will error by default, because it needs a working array size of 7.4TB.
A: If you apply SVD on the covariance matrix, your principal vectors are the same as applying SVD on the data matrix. So, mathematically they are equivalent in this case.
However, in terms of complexity, it does not make much sense to apply SVD on the covariance matrix: you have constructed the covariance matrix and then you pay for SVD which is more expensive than computing eigenvectors.
The best practice is to apply SVD directly on the data matrix, to save some flops (compared to Andrew Ng's way) and to achieve numerical stability of SVD routines (compared to eigendecomposition).
A: If anyone cares about performance difference in a specific case, I compared
SVD and Eigen-analysys on a 632x632 real symmetric matrix L using C++/Eigen.
I used BDCSVD
// compute full V, U *not* needed
Eigen::BDCSVD<Eigen::MatrixXf> svd(L,Eigen::ComputeFullV);

and SelfAdjointEigenSolver
Eigen::SelfAdjointEigenSolver<Eigen::MatrixXf> eigenSolver(N);
eigenSolver.compute(L);

and ran each block of code 100 times (clang++ -O3 -DNEBUG) on a
3.2 Ghz Intel Core i5 Mac. These are the results:
 BCSSVD ..................... 10036 milliseconds
 SelfAdjointEigenSolver ..... 10338 milliseconds

So each takes about 100 msec for a single call. So a push speed wise.
