\begin{eqnarray}
C=\sum_{i=1}^N (a_1 x_1^i + \ldots + a_n x_n^i + b – y^i)^2
\end{eqnarray}
\begin{eqnarray}
\frac {\partial C} {\partial a_m} =2\sum_{i=1}^N (a_1 x_1^i + \ldots + a_n x_n^i + b – y^i) x_m^i
\end{eqnarray}
\begin{eqnarray}
\frac{\partial C}{\partial a_m}
&=&
2
\begin{bmatrix}
x_m^1 & x_m^2 & \cdots & x_m^N
\end{bmatrix}
\begin{bmatrix}
a_1 x_1^1 + \cdots + a_n x_n^1 + b – y^1 \\
a_1 x_1^2 + \cdots + a_n x_n^2 + b – y^2 \\
\vdots \\
a_1 x_1^N + \cdots + a_n x_n^N + b – y^N
\end{bmatrix}.
\end{eqnarray}
\begin{eqnarray}
\begin{bmatrix}
\frac{\partial C}{\partial b} \\[6pt]
\frac{\partial C}{\partial a_1} \\
\vdots \\
\frac{\partial C}{\partial a_n}
\end{bmatrix}
&=&
2
\begin{bmatrix}
1 & 1 & \cdots & 1 \\
x_1^1 & x_1^2 & \cdots & x_1^N \\
\vdots & \vdots & \ddots & \vdots \\
x_n^1 & x_n^2 & \cdots & x_n^N
\end{bmatrix}
\begin{bmatrix}
b + a_1 x_1^1 + \cdots + a_n x_n^1 – y^1 \\
b + a_1 x_1^2 + \cdots + a_n x_n^2 – y^2 \\
\vdots \\
b + a_1 x_1^N + \cdots + a_n x_n^N – y^N
\end{bmatrix}.
\end{eqnarray}
\begin{eqnarray}
\begin{bmatrix}
\frac{\partial C}{\partial b} \\[6pt]
\frac{\partial C}{\partial a_1} \\
\vdots \\
\frac{\partial C}{\partial a_n}
\end{bmatrix}
&=&
2
\begin{bmatrix}
1 & 1 & \cdots & 1 \\
x_1^1 & x_1^2 & \cdots & x_1^N \\
\vdots & \vdots & \ddots & \vdots \\
x_n^1 & x_n^2 & \cdots & x_n^N
\end{bmatrix}
\left(
\begin{bmatrix}
b + a_1 x_1^1 + \cdots + a_n x_n^1 \\
b + a_1 x_1^2 + \cdots + a_n x_n^2 \\
\vdots \\
b + a_1 x_1^N + \cdots + a_n x_n^N
\end{bmatrix}
–
\begin{bmatrix}
y^1 \\
y^2 \\
\vdots \\
y^N
\end{bmatrix}
\right).
\end{eqnarray}
\begin{eqnarray}
\underbrace{
\begin{bmatrix}
\frac{\partial C}{\partial b} \\[6pt]
\frac{\partial C}{\partial a_1} \\
\vdots \\
\frac{\partial C}{\partial a_n}
\end{bmatrix}
}_{\nabla C}
&=&
2
\underbrace{
\begin{bmatrix}
1 & 1 & \cdots & 1 \\
x_1^1 & x_1^2 & \cdots & x_1^N \\
\vdots & \vdots & \ddots & \vdots \\
x_n^1 & x_n^2 & \cdots & x_n^N
\end{bmatrix}
}_{\mathbb{X}^T}
\left(
\underbrace{
\begin{bmatrix}
1 & x_1^1 & \cdots & x_n^1 \\
1 & x_1^2 & \cdots & x_n^2 \\
\vdots & \vdots & \ddots & \vdots \\
1 & x_1^N & \cdots & x_n^N
\end{bmatrix}
}_{\mathbb{X}}
\underbrace{
\begin{bmatrix}
b \\
a_1 \\
\vdots \\
a_n
\end{bmatrix}
}_{\overrightarrow{a}}
–
\underbrace{
\begin{bmatrix}
y^1 \\
y^2 \\
\vdots \\
y^N
\end{bmatrix}
}_{\overrightarrow{y}}
\right).
\end{eqnarray}
\begin{eqnarray}
\nabla C = \mathbb{X}^T ( \mathbb{X} \overrightarrow{a} – \overrightarrow{y} )
\end{eqnarray}
For an extremum point, $\nabla C = 0$, which means $\mathbb{X}^T ( \mathbb{X} \overrightarrow{a} – \overrightarrow{y} )=0$, or
\begin{eqnarray}
\overrightarrow{a} = (\mathbb{X}^T \mathbb{X})^{-1} \mathbb{X}^T \overrightarrow{y}
\end{eqnarray}