Analytical solution to linear-regression coefficient estimates I'm trying to understand matrix notation, and working with vectors and matrices. 
Right now I'd like to understand how the vector of coefficient estimates $\hat{\beta}$ in multiple regression is computed.
The basic equation seems to be  
$$
\frac{d}{d\boldsymbol{\beta}} (\boldsymbol{y}-\boldsymbol{X\beta})'(\boldsymbol{y}-\boldsymbol{X\beta}) = 0 \>.
$$
Now how would I solve for a vector $\beta$ here?
Edit: Wait, I'm stuck. I'm here now and don't know how to continue:
$    \frac{d}{d{\beta}} 
\left( 
\left(\begin{smallmatrix}
 y_1 \\
 y_2 \\
 \vdots \\
 y_n
\end{smallmatrix}\right)
-
\left(\begin{smallmatrix}
 1 & x_{11} & x_{12} & \dots & x_{1p} \\
 1 & x_{21} & x_{22} & \dots & x_{2p} \\
 \vdots & & & & \vdots \\
 1 & x_{n1} & x_{n2} & \dots & x_{np} \\
\end{smallmatrix}\right)
\left(\begin{smallmatrix}
 \beta_0 \\
 \beta_1 \\
 \vdots \\
 \beta_p
\end{smallmatrix}\right)
\right) '
\left( 
\left(\begin{smallmatrix}
 y_1 \\
 y_2 \\
 \vdots \\
 y_n
\end{smallmatrix}\right)
-
\left(\begin{smallmatrix}
 1 & x_{11} & x_{12} & \dots & x_{1p} \\
 1 & x_{21} & x_{22} & \dots & x_{2p} \\
 \vdots &  & & & \vdots \\
 1 & x_{n1} & x_{n2} & \dots & x_{np} \\
\end{smallmatrix}\right)
\left(\begin{smallmatrix}
 \beta_0 \\
 \beta_1 \\
 \vdots \\
 \beta_p
\end{smallmatrix}\right)
\right)
$
$
\frac{d}{d{\beta}} \sum_{i=1}^n \left( y_i - \begin{pmatrix} 1 & x_{i1} & x_{i2} & \dots & x_{ip} \end{pmatrix} 
\begin{pmatrix}
 \beta_0 \\
 \beta_1 \\
 \vdots \\
 \beta_p
\end{pmatrix} \right)^2$
With $x_{i0} = 1$ for all $i$ being the intercept:
$
\frac{d}{d{\beta}} \sum_{i=1}^n \left( y_i - 
\sum_{k=0}^p x_{ik} \beta_k
 \right)^2
$
Can you point me in the right direction?
 A: Here is a technique for minimizing the sum of squares in regression that actually has applications to more general settings and which I find useful. 
Let's try to avoid vector-matrix calculus altogether.
Suppose we are interested in minimizing 
$$
\newcommand{\err}{\mathcal{E}}\newcommand{\my}{\mathbf{y}}\newcommand{\mX}{\mathbf{X}}\newcommand{\bhat}{\hat{\beta}}\newcommand{\reals}{\mathbb{R}}
\err = (\my - \mX \beta)^T (\my - \mX \beta) = \|\my - \mX \beta\|_2^2 \> ,
$$
where $\my \in \reals^n$, $\mX \in \reals^{n\times p}$ and $\beta \in \reals^p$. We assume for simplicity that $p \leq n$ and $\mathrm{rank}(\mX) = p$.
For any $\bhat \in \reals^p$, we get
$$
\err = \|\my - \mX \bhat + \mX \bhat - \mX \beta\|_2^2 = \|\my - \mX \bhat\|_2^2 + \|\mX(\beta-\bhat)\|_2^2 - 2(\beta - \bhat)^T \mX^T (\my - \mX \bhat) \>.
$$
If we can choose (find!) a vector $\bhat$ such that the last term on the right-hand side is zero for every $\beta$, then we would be done, since that would imply that $\min_\beta \err \geq \|\my - \mX \bhat\|_2^2$.
But, $(\beta - \bhat)^T \mX^T (\my - \mX \bhat) = 0$ for all $\beta$ if and only if $\mX^T (\my - \mX \bhat) = 0$ and this last equation is true if and only if $\mX^T \mX \bhat = \mX^T \my$. So $\err$ is minimized by taking $\bhat = (\mX^T \mX)^{-1} \mX^T \my$.

While this may seem like a "trick" to avoid calculus, it actually has wider application and there is some interesting geometry at play.
One example where this technique makes a derivation much simpler than any matrix-vector calculus approach is when we generalize to the matrix case. Let $\newcommand{\mY}{\mathbf{Y}}\newcommand{\mB}{\mathbf{B}}\mY \in \reals^{n \times p}$, $\mX \in \reals^{n \times q}$ and $\mB \in \reals^{q \times p}$. Suppose we wish to minimize
$$
\err = \mathrm{tr}( (\mY - \mX \mB) \Sigma^{-1} (\mY - \mX \mB)^T )
$$
over the entire matrix $\mB$ of parameters. Here $\Sigma$ is a covariance matrix.
An entirely analogous approach to the above quickly establishes that the minimum of $\err$ is attained by taking
$$
\hat{\mB} = (\mX^T \mX)^{-1} \mX^T \mY \>.
$$
That is, in a regression setting where the response is a vector with covariance $\Sigma$ and the observations are independent, then the OLS estimate is attained by doing $p$ separate linear regressions on the components of the response.
A: One way which may help you understand is to not use matrix algebra, and differentiate with each respect to each component, and then "store" the results in a column vector.  So we have:
$$\frac{\partial}{\partial \beta_{k}}\sum_{i=1}^{N}\left(Y_{i}-\sum_{j=1}^{p}X_{ij}\beta_{j}\right)^{2}=0$$
Now you have $p$ of these equations, one for each beta.  This is a simple application of the chain rule:
$$\sum_{i=1}^{N}2\left(Y_{i}-\sum_{j=1}^{p}X_{ij}\beta_{j}\right)^{1}\left(\frac{\partial}{\partial \beta_{k}}\left[Y_{i}-\sum_{j=1}^{p}X_{ij}\beta_{j}\right]\right)=0$$
$$-2\sum_{i=1}^{N}X_{ik}\left(Y_{i}-\sum_{j=1}^{p}X_{ij}\beta_{j}\right)=0$$
Now we can re-write the sum inside the bracket as $\sum_{j=1}^{p}X_{ij}\beta_{j}=\bf{x}_{i}^{T}\boldsymbol{\beta}$  So you get:
$$\sum_{i=1}^{N}X_{ik}Y_{i}-\sum_{i=1}^{N}X_{ik}\bf{x}_{i}^{T}\boldsymbol{\beta}=0$$
Now we have $p$ of these equations, and we will "stack them" in a column vector.  Notice how $X_{ik}$ is the only term which depends on $k$, so we can stack this into the vector $\bf{x}_{i}$ and we get:
$$\sum_{i=1}^{N}\bf{x}_{i}\rm{Y}_{i}=\sum_{i=1}^{N}\bf{x}_{i}\bf{x}_{i}^{T}\boldsymbol{\beta}$$
Now we can take the beta outside the sum (but must stay on RHS of sum), and then take the invervse:
$$\left(\sum_{i=1}^{N}\bf{x}_{i}\bf{x}_{i}^{T}\right)^{-1}\sum_{i=1}^{N}\bf{x}_{i}\rm{Y}_{i}=\boldsymbol{\beta}$$
A: We have
$\frac{d}{d\beta} (y - X \beta)' (y - X\beta) = -2 X' (y - X \beta)$.
It can be shown by writing the equation explicitly with components. For example, write $(\beta_{1}, \ldots, \beta_{p})'$ instead of $\beta$. Then take derivatives with respect to $\beta_{1}$, $\beta_{2}$, ..., $\beta_{p}$ and stack everything to get the answer. For a quick and easy illustration, you can start with $p = 2$.
With experience one develops general rules, some of which are given, e.g., in that document.
Edit to guide for the added part of the question
With $p = 2$, we have
$(y - X \beta)'(y - X \beta) = (y_1 - x_{11} \beta_1 - x_{12} \beta_2)^2 + 
(y_2 - x_{21}\beta_1 - x_{22} \beta_2)^2$
The derivative with respect to $\beta_1$ is
$-2x_{11}(y_1 - x_{11} \beta_1 - x_{12} \beta_2)-2x_{21}(y_2 - x_{21}\beta_1 - x_{22} \beta_2)$
Similarly, the derivative with respect to $\beta_2$ is
$-2x_{12}(y_1 - x_{11} \beta_1 - x_{12} \beta_2)-2x_{22}(y_2 - x_{21}\beta_1 - x_{22} \beta_2)$
Hence, the derivative with respect to $\beta = (\beta_1, \beta_2)'$ is
$
\left( 
\begin{array}{c}
-2x_{11}(y_1 - x_{11} \beta_1 - x_{12} \beta_2)-2x_{21}(y_2 - x_{21}\beta_1 - x_{22} \beta_2) \\
-2x_{12}(y_1 - x_{11} \beta_1 - x_{12} \beta_2)-2x_{22}(y_2 - x_{21}\beta_1 - x_{22} \beta_2)
\end{array}
\right)
$
Now, observe you can rewrite the last expression as 
$-2\left( 
\begin{array}{cc}
x_{11} & x_{21} \\
x_{12} & x_{22}
\end{array}
\right)\left( 
\begin{array}{c}
y_{1} - x_{11}\beta_{1} - x_{12}\beta_2 \\
y_{2} - x_{21}\beta_{1} - x_{22}\beta_2
\end{array}
\right) = -2 X' (y - X \beta)$
Of course, everything is done in the same way for a larger $p$. 
A: You can also use formulas from Matrix cookbook. We have
$$(y-X\beta)'(y-X\beta)=y'y-\beta'X'y-y'X\beta+\beta'X'X\beta$$
Now take derivatives of each term. You might want to notice that $\beta'X'y=y'X\beta$. The derivative of term $y'y$ with respect to $\beta$ is zero. The remaining term
$$\beta'X'X\beta-2y'X\beta$$
is of form of function
$$f(x)=x'Ax+b'x,$$
in formula (88) in the book in page 11,  with $x=\beta$, $A=X'X$ and $b=-2X'y$. The derivative is given in the formula (89):
$$\frac{\partial f}{\partial x}=(A+A')x+b$$
so
$$\frac{\partial}{\partial \beta}(y-X\beta)'(y-X\beta)=(X'X+(X'X)')\beta-2X'y$$
Now since $(X'X)'=X'X$ we get the desired solution:
$$X'X\beta=X'y$$
