Finding the unbiased variance estimator in high dimensional spaces The problem comes from linear regression. Assume the regression function is linear, i.e. 
$$
f(X) = \beta_0+\sum_{j=1}^pX_j\beta_j
$$
.Given a set of training data $(x_1, y_1),\ldots,(x_N,y_N)$,we try to estimate the parameters $\beta$ by minimizing the residual sum of squares: 
$$\text{RSS}(\beta) = (\sum_{i=1}^N(y_i-\beta_0-\sum_{j=1}^px_{ij}\beta_j)^2$$
Here each $x_i=(x_{i1},\ldots,x_{ip})^T$ is vector in $\mathbb{R}^p$,$y_i\in\mathbb{R}$,and $\beta=(\beta_0,\ldots,\beta_{p+1})$. Let $\mathbf{X}$ denote matrix $(x_1,\ldots,x_N)^T$, $\mathbf{y}=(y_1,\ldots,y_N)$ and assume it has full column rank, it's easy to get
$$\hat{\beta}=(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}$$
Suppose $\mathbf{\hat{y}}=\mathbf{X}\hat{\beta}$, and $\mathbf{y}=\mathbf{X}\beta$ for some $\beta$, to pin down the sampling properties of $\hat{\beta}$, now assume that the observations $y_i$ are uncorrelated and have constant variance $\sigma^2$, and that the $x_i$ are fixed(non random). It's easy to see that
$$\text{Var}(\hat\beta) = (\mathbf{X}^T\mathbf{X})^{-1}\sigma^2$$
,so we turn to estimating the $\sigma^2$, and now comes the question:
Let $\hat\sigma$ defined as below:
$$\hat\sigma^2=\frac{1}{K}\sum_{i=1}^N(y_i-\hat{y_i})^2$$
It is said that only when $K=N-p-1$ did $\hat\sigma$ be an unbiased estimator of $\sigma$. But how do $p$(the dimension of $x$) be introduced in?
 A: Well it might be an overkill but I think this proof is OK. I would use basic linear algebra tools. Starting with slight change in your notations:
Let $X$ denote the matrix $(x_1\; x_2 \;...\; x_N)^T$ where $x_i=(1 \; x_{i2} \; x_{i3} \; ... \; x_{ip})$. So now we have $p-1$ covariates. Our model is 
$$
y=X\beta+\epsilon
$$
with $E(\epsilon)=0$ and $Var(\epsilon)=\sigma^2I$
Now, we can write 
$$
\hat{y}=X\hat{\beta}=P_Xy
$$
where $P_X$ is the projection matrix onto the space spanned by the columns of $X$. That is,
$$
P_X=X(X^TX)^{-1}X^T.
$$
Note that $P_X$ is $n\times n$ matrix.
Now, note that we can actually write $\sum_{i=1}^{N}(y_i-\hat{y}_i)^2$ as $||y-P_Xy||^2 $ with $||\cdot||$ being the $\ell_2$ norm of a vector.
Next,
$$
||y-P_Xy||^2=||(I-P_X)y||^2=y^T(I-P_X)y
$$
where the last equality holds since $I-P_X$ is a projection matrix.
Note also that $(I-P_X)X=0$ and hence $y^T(I-P_X)y=\epsilon^T(I-P_X)\epsilon$ (using our model definition).
If we assume normality of the error, we could continue to derive the appropriate $\chi^2$ distribution. 
But you wrote nothing regarding the distribution so I am now turning to the expectation
Recall that $I-P_X$ is $N\times N$ projection matrix. We can use the eigendecomposition of a matrix to write
$$
(I-P_X)=U\Lambda U^T
$$
where $UU^T=I$ and $\Lambda$ is a diagonal matrix with the diagonal being the eigenvalues of $I-P_X$. 
Since $I-P_X$ is a projection matrix, all of its eigenvalues are equal to zero or one. How many "ones" there are? as the dimension of the subspace that $I-P_X$ projecting onto. Since this is the complementary of the subspace spanned by $P_X$, we get $n-p$ ones in the diagonal, and by using the fact that $UU^T=I$ we get the final answer
which is 
$$
E(||(I-P_X)y||^2)= E(\epsilon^TU\Lambda U^T\epsilon)=\sum_{j=1}^n\lambda_j\sigma^2=(n-p)\sigma^2
$$
and as written before there are $n-p \; \lambda_j$'s that are equal to one and $p \; \lambda_j$'s that equals to zero  ($\lambda_j, i=1,...,n$ are the eigenvalues of $I-P_X$).
Some relevant Wikipedia pages for the linear algebra:
http://en.wikipedia.org/wiki/Projection_(linear_algebra)
http://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix
http://en.wikipedia.org/wiki/Symmetric_matrix
