Why are principal components in PCA (eigenvectors of the covariance matrix) mutually orthogonal? Why are principal components in PCA mutually orthogonal? I know that PCA can be calculated by eig(cov(X)), where X is centered. But I do not see why the eigenvectors should be orthogonal.
 A: I think it might help to pull back from the mathematics and think about the goal of PCA. In my mind, PCA is used to represent large-dimensional data sets (many variables) in the clearest way possible--i.e. the way that reveals as much of the underlying data structure as possible. 
For an example, let's consider a data set with 3 variables: height, weight and age of people. If you have this data, you could create 3 separate 2-dimensional scatter plots, Height vs. Weight, Weight vs. Age, and Age vs. Height--and if you graphed them, they would all look like 'squares' with base=X-axis and height=Y-axis. 
But if you suspect there is an interesting relationship among all 3 variables, you would probably want to create a 3-dimensional plot of the data, with X, Y and Z axes--which would create a 3-D 'cube' plot. This plot might reveal an interesting relationship that you would like to communicate to people. But of course, you can't print a 3-dimensional plot, so you have to project the data onto a 2-dimensional piece of paper, which means you have to choose 2 dimensions to prioritize, at the expense of the third (which would be orthogonal to the piece of paper you are printing on). 
You could try to include the 3rd dimension visually with the use of color-coding or various sized bubbles for the points. Or, you could rotate the plot (in your mind, or with software) until you find a new projection that expresses as much of the 3-D information in 2-D space as possible--think visually of just rotating the 'cube' around until the underlying data relationship you want to show is a clear as possible for printing on paper. As you rotate the 3-D 'cube' through 2-D space, you create new synthetic axes, and these new synthetic axes are orthogonal to each other--and correspond to the length and width of the paper you are printing on. If you travel along these new synthetic axes, you are moving through multiple dimensions of the original data (height, weight and age) at the same time, but you can travel along the new synthetic X-axis (the width of the paper) without moving along the synthetic Y-axis. 
We can think of this visually, because our brains understand 3-dimensional spaces, but things quickly become problematic if you are talking about higher dimensional data sets. We can't imagine 9-dimensional 'hyper-cubes' (at least I can't), but we often have to deal with data sets that contain many variables. We can use software (or grueling math) to 'rotate' the 9-dimensional data through space, until we find the new projection that represents as much of the higher-dimensional data structure as possible--for printing on a 2-D page. 
This is exactly what PCA does. Again, for simplicity, consider the earlier 3-D data set example plotted in a 'cube' space--we would see something like a cloud of points. PCA simply rotates that cloud until it finds the 'longest' straight line possible through that cloud--the direction of this line becomes PC1. Then, with PC1 fixed, the data point cloud is rotated again, along PC1, until the next 'longest' orthogonal axis is found, which is PC2. Then you can print a new 2-D plot (PC1 vs. PC2) that captures as much of the 3-D data structure as possible. And of course you can keep going and find PC3, PC4 and so on, if it helps understand the data.
Then, the PCA results will tell you how much of the data variance is explained by the new synthetic principal components, and if the PCA axes capture more of the data variance than you would expect to occur by random chance, we can infer that there is a meaningful relationship among the original measured variables.
All the discussion about eigenvectors and matrix algebra is a little bit beside the point in my opinion (and also, I'm not that mathematically inclined)--orthogonal axes are just an inherent part of this type of matrix algebra. So, citing the mathematical foundations of orthogonal axes doesn't really explain why we use this approach for PCA. We use this matrix algebra in statistical analysis because it helps us reveal important characteristics of data structures that we are interested in.
A: The covariance matrix is symmetric. If a matrix $A$ is symmetric, and has two eigenvectors $u$ and $v$, consider $Au = \lambda u$ and $Av = \mu v$.
Then by symmetry (and writing $'$ for transpose):
$$u'Av = u'A'v = (Au)'v = \lambda u'v$$
More directly:
$$u'Av = u'(\mu v) = \mu u'v$$
Since these are equal we obtain $(\lambda - \mu)u'v =0$. So either $u'v = 0$ and the two vectors are orthogonal, or $\lambda - \mu = 0$ and the two eigenvalues are equal. In the latter case, the eigenspace for that repeated eigenvalue can contain eigenvectors which are not orthogonal. So your instinct to question why the eigenvectors have to be orthogonal was a good one; if there are repeated eigenvalues they may not be! What if your sample covariance is the identity matrix? This has repeated eigenvalue $1$ and any two non-zero vectors are eigenvectors, orthogonal or not. (Thinking out such special cases is often a good way to spot counter-examples.)
If a symmetric matrix has a repeated eigenvalue, we can choose to pick out orthogonal eigenvectors from its eigenspace. That's what we want to do in PCA, because finding orthogonal components is the whole point of the exercise. Of course it's unlikely that your sample covariance matrix will have repeated eigenvalues - if so, it would only have taken a small perturbation of your data to make them unequal - but we should take care to define our algorithm so it really does pick out orthogonal eigenvectors. (Note that which it picks out, and in what order, is arbitrary. Think back to the identity matrix and all its possible orthogonal sets of eigenvectors! This is a generalisation of the arbitrary choice between $v$ and $-v$ for a unique eigenvalue. Output from two different implementations of PCA may look quite different.)
To see why it's guaranteed that we can set up our implementation of eig(cov(X)) in this manner - in other words, why there are always "just enough" orthogonal vectors we can pick out from that eigenspace - you need to understand why geometric and algebraic multiplicities are equal for symmetric matrices. If the eigenvalue appears twice we can pick out two orthogonal eigenvectors; thrice and we can pick out three, and so on. Several approaches are raised in this mathematics stack exchange thread but the usual method is via the Schur decomposition. The result you are after is probably proved in your linear algebra textbook as the "spectral theorem" (though that phrase can also refer to several more general results) or perhaps under a more specific name like "symmetric eigenvalue decomposition". Symmetric matrices have several nice properties that it's worth knowing, e.g. their eigenvalues are real, so we can find real eigenvectors, with obvious implications for PCA.
Finally, how can we write an implementation that achieves this? I will consider two implementations of PCA in R. We can see the code for princomp: look at methods(princomp) then getAnywhere(princomp.default) and we observe edc <- eigen(cv, symmetric = TRUE). So eigen will use LAPACK routines for symmetric matrices. Checking the LAPACK Users' Guide (3rd edition) for "symmetric eigenproblems" we see it firstly decomposes $A = QTQ'$ where $T$ is symmetric tridiagonal and $Q$ is orthogonal, then decomposes $T = S \Lambda S'$ where $\Lambda$ is diagonal and $S$ orthogonal. Then writing $Z = QS$ we have diagonalized $A = Z \Lambda Z'$. Here $\Lambda$ is the vector of eigenvalues (of $T$ and also of $A$ - they work out the same) and since $Z$ is the product of two orthogonal matrices it is also orthogonal. The computed eigenvectors are the columns of $Z$ so we can see LAPACK guarantees they will be orthonormal (if you want to know quite how the orthogonal vectors of $T$ are picked, using a Relatively Robust Representations procedure, have a look at the documentation for DSYEVR). So that's one approach, but for numerical reasons it'd be better to do a singular value decomposition. If you look under the bonnet of another PCA function in R, you'll see this is how prcomp works. R uses a different bunch of LAPACK routines to solve this problem.
