I can give you my own explanation/proof of the PCA, which I think is really simple and elegant, and doesn't require anything except basic knowledge of linear algebra. It came out pretty lengthy, because I wanted to write in simple accessible language.
Suppose we have some $M$ samples of data from an $n$-dimensional space. Now we want to project this data on a few lines in the $n$-dimensional space, in a way that retains as much variance as possible (that means, the variance of the projected data should be as big compared to the variance of original data as possible).
Now, let's observe that if we translate (move) all the points by some vector $\beta$, the variance will remain the same, since moving all points by $\beta$ will move their arithmetic mean by $\beta$ as well, and variance is lineary proportional to $\sum_{i=1}^M ||x_i - \mu||_2$. Hence we translate all the points by $-\mu$, so that their arithmetic mean is $0$, for computational comfort. Let's denote the translated points as $x_i' = x_i - \mu$. Let's also observe, that the variance can be now expressed simply as $\sum_{i=1}^M ||x_i'||_2$.
Now the choice of the line. We can describe any line as set of points that satisfy the equation $x = \alpha * v + w$, for some vectors $v,w$. Note that if we move the line by some vector $\gamma$ orthogonal to $v$, then all the projections on the line will also be moved by $\gamma$, hence the mean of the projections will be moved by $\gamma$, hence the variance of the projections will remain unchanged. That means we can move the line parallel to itself, and not change the variance of projections on this line. Again for convenience purposes let's limit ourselves to only the lines passing through the zero points (this means lines described by $x = \alpha *v$).
Alright, now suppose we have a vector $v$ that describes the direction of a line that is a possible candidate for the line we search for. We need to calculate variance of the projections on the line $\alpha * v$. What we will need are projection points and their mean. From linear algebra we know that in this simple case the projection of $x_i'$ on $\alpha * v$ is $\frac{<x_i, v>}{||v||_2}$. Let's from now on limit ourselves to only unit vectors v. That means we can write the length of projection of point $x_i'$ on $v$ simply as $<x_i', v>$.
In some of the previous answers someone said that PCA minimizes the sum of squares of distances from the chosen line. We can now see it's true, because sum of squares of projections plus sum of squares of distances from the chosen line is equal to the sum of squares of distances from point $0$. By maximizing the sum of squares of projections, we minimize the sum of squares of distances and vice versa, but this was just a thoughtful digression, back to the proof now.
As for the mean of the projections, let's observe that $v$ is part of some orthogonal base of our space, and that if we project our data points on every vector of that base, their sum will cancel out (it's like that because projecting on the vectors from the base is like writing the data points in the new orthogonal base). So the sum of all the projections on the vector $v$ (let's call the sum $S_v$) and the sum of projections on other points from the base (let's call it $S_o$) is 0, because it's the mean of the data points. But $S_v$ is orthogonal to $S_o$! That means $S_o = S_v = 0$.
So the mean of our projections is $0$? Well, that's convenient, because that means the variance is just the sum of squares of lengths of projections, or in symbols $\sum_{i=1}^M (x_i' \cdot v)^2 = \sum_{i=1}^M v^T \cdot x_i'^T \cdot x_i' \cdot v = v^T \cdot (\sum_{i=1}^M x_i'^T \cdot x_i) \cdot v$.
Well well, suddenly the covariance matrix popped out. Let's denote it simply by $X$. It means we are now looking for a unit vector $v$, that maximizes $v^T \cdot X \cdot v$, for a semi-positive definite matrix $X$.
Now, let's take the eigenvectors and eigenvalues of matrix $X$, and denote them by $e_1, e_2, \dots , e_n$ and $\lambda_1 , \dots, \lambda_n$ respectively, such that $\lambda_1 \geq \lambda_2 , \geq \lambda_3 \dots $. If the values $\lambda$ do not duplicate, eigenvectors form an orthonormal base. If they do, we choose the eigenvectors in a way that they form an orthonormal base.
Now let's calculate $v^T \cdot X \cdot v$ for eigenvector $e_i$. We have $e_i^T \cdot X \cdot e_i = e_i^T \cdot (\lambda_i * e_i) = \lambda_i (||e_i||_2)^2 = \lambda_i$.
Pretty good, this gives us $\lambda_1$ for $e_1$. Now let's take an arbitrary vector $v$. Since eigenvectors form an orthonormal base, we can write $v = \sum_{i=1}^n e_i * <v, e_i>$, and we have $\sum_{i=1}^n <v, e_i>^2 = 1$. Let's denote $\beta_i = <v, e_i>$.
Now let's count $v^T \cdot X \cdot v$. We rewrite $v$ as a linear combination of $e_i$, and get : $(\sum_{i=1}^n \beta_i * e_i)^T \cdot X \cdot (\sum_{i=1}^n \beta_i *e_i) = (\sum_{i=1}^n \beta_i * e_i) \cdot (\sum_{i=1}^n \lambda_i * \beta_i * e_i) = \sum_{i=1}^n \lambda_i *(\beta_i)^2 * (||e_i||_2)^2$.
The last equation comes from the fact thet eigenvectors where chosen to be pairwaise orthogonal, so their dot product's are zero. Now, because all eigenvectors are also unit, we can write $v^T \cdot X \cdot v = \sum_{i=1}^n \lambda_i * \beta_i^2$, where $\beta_i ^2$ are all positive, and sum to $1$.
That means, that the variance of the projections is a weighted mean of eigenvalues. Certainly, it is always less then the biggest eigenvalue, which is why it should be our choice of the first vector.
Now suppose we want another vector. We should chose it from a space orthogonal to the already chosen one, that means the subspace $lin(e_2, e_3, \dots , e_n)$. By analogical inference we arrive at the conclusion, that the best vector to project on is $e_2$. And so on, and so on...
Btw. it should be now clear, why the variance retained can be expresed by $\frac{\sum_{i=1}^k \lambda_i}{\sum_{i=1}^n \lambda_i}$.
Hope this helps and I didn't make any big mistakes.
We should also justify the greedy choice of vectors. When we want too choose $k$ vectors to project onto, it might not be the best idea to first choose the best vector, then the best from what reamains and so on. I'd like to argue that in this case it is justified and makes no difference. Lets denote the $k$ vector we wish to project onto by $v_1, \dots , v_k$. Also, let's assume the vectors are pairwaise orthogonal. As we already know, the total variance of the projections on those vectors can be expressed by : $\sum_{j=1}^k \sum_{i=1}^n \lambda_i * \beta_{ij}^2 = \sum_{i=1}^n \lambda_i * \gamma_i$ where $\gamma_i = \sum_{j=1}^k \beta_{ij}^2$.
Now, let's write e_i in some orthonormal base that includes $v_1, \dots , v_k$. Let's denote the rest of the base as $u_1, \dots, u_{n-k}$. We can see that $e_i = \sum_{j=1}^k \beta_{ij} v_j + \sum_{j=1}^{n-k} \theta_j * <e_i, u_j>$.
Because $||e_i||_2 = 1$, we have $\sum_{j=1}^k \beta_{ij}^2 + \sum_{j=1}^{n-k} \theta_j^2 = 1$, and hence $\gamma_i \leq 1$ for all $i$.
Now we have a similar case to one vector only, we now know that the total variance of projections is $\sum_{i=1}^n \lambda_i * \gamma_i$ with $\gamma_i \leq 1$ and $\sum_{i=1}^n \gamma_i = k$. This is yet another weighted mean, and is certainly no more than $\sum_{i=1}^k \lambda_i$ which corresponds to projecting on $k$ eigenvectors corresponding to biggest eigenvalues.
