Given a PCA (or SVD) approximation of matrix $X$ with a matrix $\hat X$, we know that $\hat X$ is the best low-rank approximation of $X$.
Is this according to the induced $\parallel \cdot \parallel_2$ norm (i.e. the largest eigenvalue norm) or according to the Frobenius $\parallel \cdot \parallel_F$ norm ?
Let's start with defining the norms. For a matrix $X$, operator $2$-norm is defined as $$\|X\|_2 = \mathrm{sup}\frac{\|Xv\|_2}{\|v\|_2} = \mathrm{max}(s_i)$$ and Frobenius norm as $$\|X\|_F = \sqrt {\sum_{ij} X_{ij}^2} = \mathrm{tr}(X^\top X) = \sqrt{\sum s_i^2},$$ where $s_i$ are singular values of $X$, i.e. diagonal elements of $S$ in the singular value decomposition $X = USV^\top$.
PCA is given by the same singular value decomposition when the data are centered. $US$ are principal components, $V$ are principal axes, i.e. eigenvectors of the covariance matrix, and the reconstruction of $X$ with only the $k$ principal components corresponding to the $k$ largest singular values is given by $X_k = U_k S_k V_k^\top$.
The Eckart-Young theorem says that $X_k$ is the matrix minimizing the norm of the reconstruction error $\|X-A\|$ among all matrices $A$ of rank $k$. This is true for both, Frobenius norm and the operator $2$-norm. As pointed out by @cardinal in the comments, it was first proved by Schmidt (of Gram-Schmidt fame) in 1907 for the Frobenius case. It was later rediscovered by Eckart and Young in 1936 and is now mostly associated with their names. Mirsky generalized the theorem in 1958 to all norms that are invariant under unitary transformations, and this includes the operator 2-norm.
This theorem is sometimes called Eckart-Young-Mirsky theorem. Stewart (1993) calls it Schmidt approximation theorem. I have even seen it called Schmidt-Eckart-Young-Mirsky theorem.
Let $X$ be of full rank $n$. As $A$ is of rank $k$, its null space has $n-k$ dimensions. The space spanned by the $k+1$ right singular vectors of $X$ corresponding to the largest singular values has $k+1$ dimensions. So these two spaces must intersect. Let $w$ be a unit vector from the intersection. Then we get: $$\|X-A\|^2_2 \ge \|(X-A)w\|^2_2 = \|Xw\|^2_2 = \sum_{i=1}^{k+1}s_i^2(v_i^\top w)^2 \ge s_{k+1}^2 = \|X-X_k\|_2^2,$$ QED.
We want to find matrix $A$ of rank $k$ that minimizes $\|X-A\|^2_F$. We can factorize $A=BW^\top$, where $W$ has $k$ orthonormal columns. Minimizing $\|X-BW^\top\|^2$ for fixed $W$ is a regression problem with solution $B=XW$. Plugging it in, we see that we now need to minimize $$\|X-XWW^\top\|^2=\|X\|^2-\|XWW^\top\|^2=\mathrm{const}-\mathrm{tr}(WW^\top X^\top XWW^\top)\\=\mathrm{const}-\mathrm{const}\cdot\mathrm{tr}(W^\top\Sigma W),$$ where $\Sigma$ is the covariance matrix of $X$, i.e. $\Sigma=X^\top X/(n-1)$. This means that reconstruction error is minimized by taking as columns of $W$ some $k$ orthonormal vectors maximizing the total variance of the projection.
It is well-known that these are first $k$ eigenvectors of the covariance matrix. Indeed, if $X=USV^\top$, then $\Sigma=VS^2V^\top/(n-1)=V\Lambda V^\top$. Writing $R=V^\top W$ which also has orthonormal columns, we get $$\mathrm{tr}(W^\top\Sigma W)=\mathrm{tr}(R^\top\Lambda R)=\sum_i \lambda_i \sum_j R_{ij}^2 \le \sum_{i=1}^k \lambda_k,$$ with maximum achieved when $W=V_k$. The theorem then follows immediately.
See the following three related threads:
This proof I found somewhere online but it is wrong (contains a gap), as explained by @cardinal in the comments.
Frobenius norm is invariant under unitary transformations, because they do not change the singular values. So we get: $$\|X-A\|_F=\|USV^\top - A\| = \|S - U^\top A V\| = \|S-B\|,$$ where $B=U^\top A V$. Continuing: $$\|X-A\|_F = \sum_{ij}(S_{ij}-B_{ij})^2 = \sum_i (s_i-B_{ii})^2 + \sum_{i\ne j}B_{ij}^2.$$ This is minimized when all off-diagonal elements of $B$ are zero and all $k$ diagonal terms cancel out the $k$ largest singular values $s_i$ [gap here: this is not obvious], i.e. $B_\mathrm{optimal}=S_k$ and hence $A_\mathrm{optimal} = U_k S_k V_k^\top$.
