I'll post this proof of least squares as this seems appropriate here. It's least squares from a linear algebra point of view, and adapted from Friedberg's Linear Algebra. It requires some more involved linear algebra arguments, but I think it gives a good perspective on least squares. Plus, the normal equations just fall right out of this derivation.
First, we assume we have $N$ data points $(y_1, \mathbf{x}_1), ..., (y_N, \mathbf{x}_N),$ where $\mathbf{x}_i = (1, x_{i1}, ..., x_{ip}).$ Here $y_i$ and $x_{ij}$ are taken to be real numbers. Letting $\mathbf{X}$ be the matrix whose $i^\text{th}$ row is $\mathbf{x}_i$ and $\bf{y}$ the vector whose $i^\text{th}$ component is $y_i$, the problem of least squares can be stated as finding $\beta$ satisfying:
$$
\underset{\beta \in \mathbb{R}^{p + 1}}{\text{argmin}} \ ||\mathbf{y} - \mathbf{X}\beta||^2.
$$
Remember from linear algebra, that $\bf{X}$ represents a linear transformation $\mathbb{R}^{p + 1} \to \mathbb{R}^N$ given by $\mathbf{w} \mapsto \mathbf{X}\mathbf{w},$ so it makes sense to define $W$ to be the range of $\bf{X}.$ Now, from linear algebra we know that we can write
$$
\mathbf{y} = \mathbf{u} + \mathbf{z}
$$
for some $\mathbf{u} \in W, \mathbf{z} \in W^\perp,$ where $W^\perp$ is the subspace of $\mathbb{R}^{N}$ orthogonal to $W$, or the orthogonal complement of $W.$ This vector $\mathbf{u}$ is what is usually called the orthogonal projection of $\bf{y}$ onto $W$, and already satisfies:
$$
||\mathbf{y} - \mathbf{u}||^2 \leq ||\mathbf{y} - \mathbf{w}||^2
$$
for all $\mathbf{w} \in W.$ Now, since $\mathbf{u} \in W$, then by definition $\mathbf{u} = \mathbf{X}\hat{\beta}$ for some $\hat{\beta} \in \mathbb{R}^{p + 1}$ (since $W$ is the range of $\mathbf{X}$), and
$$
\mathbf{X}\hat{\beta} - \mathbf{y} = -\mathbf{z} \in W^\perp
$$
since $W^\perp$ is a subspace. Since $W^\perp$ is the orthogonal complement of $W,$ then for any $\alpha \in \mathbb{R}^{p + 1}$ (again since $W$ is the range of $\bf{X}$):
$$
0 = \left<\mathbf{X}\hat{\beta} - \mathbf{y}, \mathbf{X}\alpha\right> = \left<\mathbf{X}^*(\mathbf{X}\hat{\beta} - \mathbf{y}), \alpha\right>.
$$
Here $\mathbf{X}^*$ is called the adjoint of $\mathbf{X},$ and since $\mathbf{X}$ only has real entries, $\mathbf{X}^* = \mathbf{X}^T.$ The last equation implies that
$$
0 = \mathbf{X}^*(\mathbf{X}\hat{\beta} - \mathbf{y}) = \mathbf{X}^T(\mathbf{X}\hat{\beta} - \mathbf{y}) = \mathbf{X}^T\mathbf{X}\hat{\beta} - \mathbf{X}^T\mathbf{y}.
$$
If you assume $\mathbf{X}^T\mathbf{X}$ is invertible, then this gives you:
$$
\hat{\beta} = \left(\mathbf{X}^T\mathbf{X}\right)^{-1}\mathbf{X}^T\mathbf{y}
$$
which are precisely the normal equations! You see that arguing in this manner bypasses the messy summation notation, and also generalizes least squares to the case where $\mathbf{y}_i, \mathbf{x}_i$ (it is easy to extend this to the case where you are regressing on vectors rather than just numeric values) are both vectors of arbitrary dimension (not necessarily equal dimension).