PCA eigenvectors I came across many different derivations of PCA, but still I have some difficulties with fully understanding this algorithm. So, assume we have covariance symmetric matrix of our data $A\in \mathbb R^{n\times n}$ and our PCA task is to maximize variance of projections onto some vector $u$  with constraint $u^T u=1$, which is given by(after some derivations): $$f(u)=u^T Au$$ After using Largrange maximization, we know that  expression below maximize above function: $$Au=\lambda u$$
and as a result optimal value of $f(u)$ is given by: $$f(u)=u^T u \lambda = \lambda$$ so our function $f(u)$ has maximal value for eigenvector $u$ with largest eigenvalue $\lambda$. And then we can choose  for example two largest eigenvalues in order to have projections of our data onto two dimensional space which is spanned by corresponding orthogonal eigenvectors which explains huge part of variance in our data due to $\operatorname{tr}A=\sum \lambda_i$
My question is following: Why we can choose arbitrarily number of eigenvectors and form them as a basis onto which we project our original data?  Why those eigenvectors with largest variance given by corresponding eigenvalues form a basis(our variance maximizing function $f(u)$ is determined by vector $u$ not a matrix $U$)?  Derivations of PCA dont seek to maximize variance of projections onto matrix $U$ but only onto vector $u$. So generally how we generalize our formula for maximization of projections onto vector $u(f(u)$ to maximization of variance of projections onto basis given by $U(f(U))$ which is formed from orthogonal eigenvectors? How you can explain  that in every PCA derivation, we maximize function of one particular vector ($f(u)$) and finally we can choose arbitrarily number of $n$ eigenvectors with largest eigenvlue to project our data onto $\mathbb R^n$ space.
Sorry for few questions instead of one, but I try to express my concerns best as possible. 
 A: We may derive PCA by finding directions that maximize the variance, and eigen vectors naturally arise as solutions. Than we show that any linear combination of these solutions is itself a solution, thus forming a basis spanning the solution space.
$$\operatorname{Var}_u(x)=\operatorname{E}((x\cdot u-\operatorname{E}(x\cdot u))^2)$$
$$\operatorname{E}((u(\cdot (x-\operatorname{E}(x\cdot u)))^2)$$
$$u\cdot \operatorname{E}((x-\operatorname{E}(x))(x-\operatorname{E}(x))^t)\cdot u$$
$$\implies$$
$$\operatorname{Var}_u(x)=u^tCu$$
Let
$$v=U\cdot c$$
and
$$u=\sum_{i=0}^n c_i u_i$$
be the projection of a vector $v$ onto the subspace spanned by the first $n$ basis vectors $U=\{u_i\}$. 
This addresses your question about projections onto a subspace, we are computing the entire basis $U$ at once. We employ Lagrange multipliers to constrain our solutions $u_i$ to be of unit length.
$$L=\left(\sum_{i=0}^nc_iu_i\right)^tC\left(\sum_{i=0}^nc_iu_i\right)+\sum_{i=0}^n\lambda_i(1-u_i^tu_i)$$
We require all partial derivatives to be zero.
$$\frac{\partial}{\partial u_i}L=(C^t+C)c^2_iu_i-\lambda_i u_i$$
$$C^t=C$$
$$0=2Cc^2_iu_i-\lambda_i u_i$$
$$Cu_i=\frac{\lambda_i}{2c^2_i}u_i$$
We recognize the familiar eigen-equation, with $\frac{\lambda_i}{2c^2_i}$ being the eigenvalue corresponding to the eigen-vector $u_i$. 
Thus, a projection of a vector onto the subspace spanned by particular solutions $U=\{u_i\}$ (eigen-basis) is itself a solution, maximizing the lagrangian.
