[proofplan] We decompose the centered response vector $y-\bar{y}\mathbb{1}_n$ into the fitted centered part $\hat{y}-\bar{y}\mathbb{1}_n$ and the residual part $y-\hat{y}$. The intercept hypothesis places the fitted centered part in $\operatorname{Range}(X)$, while the normal equations show that the residual is orthogonal to $\operatorname{Range}(X)$. The Euclidean Pythagorean identity then gives the desired equality of squared norms. [/proofplan] [step:Decompose the centered response into fitted and residual parts] Because $X$ has full column rank, the matrix $X^\top X \in \mathbb{R}^{p \times p}$ is invertible, so $H = X(X^\top X)^{-1}X^\top$ is well-defined. Define \begin{align*} a &:= \hat{y} - \bar{y}\mathbb{1}_n \in \mathbb{R}^n,\\ b &:= y - \hat{y} \in \mathbb{R}^n. \end{align*} Then \begin{align*} y - \bar{y}\mathbb{1}_n = (\hat{y} - \bar{y}\mathbb{1}_n) + (y - \hat{y}) = a + b. \end{align*} Also $\hat{y} = Hy = X((X^\top X)^{-1}X^\top y)$ belongs to $\operatorname{Range}(X)$. Since $\mathbb{1}_n \in \operatorname{Range}(X)$ by hypothesis and $\operatorname{Range}(X)$ is a linear subspace of $\mathbb{R}^n$, we have \begin{align*} a = \hat{y} - \bar{y}\mathbb{1}_n \in \operatorname{Range}(X). \end{align*} [/step] [step:Show that the residual is orthogonal to the fitted centered part] Let $z \in \operatorname{Range}(X)$. Then there exists $c \in \mathbb{R}^p$ such that $z = Xc$. Using $b = y - Hy = (I_n - H)y$, where $I_n \in \mathbb{R}^{n \times n}$ is the identity matrix, we compute \begin{align*} z^\top b &= (Xc)^\top (I_n - H)y\\ &= c^\top X^\top(I_n - H)y\\ &= c^\top\left(X^\top - X^\top X(X^\top X)^{-1}X^\top\right)y\\ &= c^\top(X^\top - X^\top)y\\ &= 0. \end{align*} Thus $b$ is orthogonal to every vector in $\operatorname{Range}(X)$. Since $a \in \operatorname{Range}(X)$, it follows that \begin{align*} a^\top b = 0. \end{align*} [guided] The residual vector is $b = y-\hat{y} = y-Hy = (I_n-H)y$. To prove the ANOVA identity, we need the two summands in the decomposition \begin{align*} y-\bar{y}\mathbb{1}_n = a+b \end{align*} to be orthogonal. Since the previous step showed that $a \in \operatorname{Range}(X)$, it is enough to prove that the residual $b$ is orthogonal to every vector in $\operatorname{Range}(X)$. Let $z \in \operatorname{Range}(X)$. By definition of the range, there is a vector $c \in \mathbb{R}^p$ such that $z = Xc$. Then \begin{align*} z^\top b &= (Xc)^\top (I_n - H)y\\ &= c^\top X^\top(I_n - H)y. \end{align*} Now substitute the definition $H = X(X^\top X)^{-1}X^\top$: \begin{align*} X^\top(I_n - H) &= X^\top - X^\top X(X^\top X)^{-1}X^\top\\ &= X^\top - X^\top\\ &= 0. \end{align*} Therefore $z^\top b = 0$. Since this holds for every $z \in \operatorname{Range}(X)$ and $a \in \operatorname{Range}(X)$, we obtain $a^\top b = 0$. [/guided] [/step] [step:Apply the Euclidean Pythagorean identity to obtain the ANOVA decomposition] Using $y-\bar{y}\mathbb{1}_n = a+b$ and $a^\top b = 0$, we expand the squared Euclidean norm: \begin{align*} |y-\bar{y}\mathbb{1}_n|^2 &= |a+b|^2\\ &= (a+b)^\top(a+b)\\ &= a^\top a + 2a^\top b + b^\top b\\ &= |a|^2 + |b|^2\\ &= |\hat{y}-\bar{y}\mathbb{1}_n|^2 + |y-\hat{y}|^2. \end{align*} By the definitions of $TSS$, $ESS$, and $RSS$, this is exactly \begin{align*} TSS = ESS + RSS. \end{align*} [/step]