$x_i$は$p-1$個の要素を持つ説明変数ベクトルとする。
$y_i$は目的変数とすると線形モデルは下記の様に表される
y_i=\beta_1+\sum_{j=2}^{p}x_{ij}\beta_j
この時、$\beta$を切片項を含むp要素のベクトルとし
$X_1$ ~ $X_N$を切片項の1を含むN行p列の行列としてまとめると
\hat{y}=X^T\hat{\beta}
とかける。この時、
RSS(\beta)=\sum_{i=1}^N(y_i -x_i^T\beta)^2
を最小とする$\beta$を求めたい。($X\beta)^Tyとy^T(X\beta)$はベクトルであることからかける順番は気にしなくてよく
\begin{align}
RSS(\beta)&=(y-X\beta)^T(y-X\beta)\\
&=y^Ty-(X\beta)^Ty-y^TX\beta+(X\beta)^TX\beta\\
&=y^Ty-2(X\beta)^Ty+(X\beta)^TX\beta\\
&=y^Ty-2(X\beta)^Ty+\beta^TX^TX\beta
\end{align}
この値を最小化する$\beta$を求めたい。$f(\beta)$とみて$\beta$について偏微分し0となる$\beta$を求める。まずは第二項の$P(\beta)=2(X\beta)^Ty$について考える
\qquad P(\beta) = 2
\begin{bmatrix}
\begin{pmatrix}
x_{11} & x_{12} & \cdots & x_{1p}\\
x_{21} & x_{22} & \cdots & x_{2p}\\
\cdots & \cdots & \cdots & \cdots\\
x_{N1} & x_{N2} & \cdots & x_{Np}
\end{pmatrix}
\begin{pmatrix}
\beta_{1}\\
\beta_{2}\\
\cdots\\
\beta_{p}\\
\end{pmatrix}
\end{bmatrix}^T
\begin{pmatrix}
y_{1}\\
y_{2}\\
\cdots\\
y_{N}\\
\end{pmatrix}\\
P(\beta)=2
\begin{bmatrix}
\begin{pmatrix}
x_{11}\beta_{1} + \cdots + x_{1p}\beta_{p}\\
x_{21}\beta_{1} + \cdots + x_{2p}\beta_{p}\\
\cdots + \cdots + \cdots\\
x_{N1}\beta_{1} + \cdots + x_{Np}\beta_{p}\\
\end{pmatrix}
\end{bmatrix}^T
\begin{pmatrix}
y_{1}\\
y_{2}\\
\cdots\\
y_{N}\\
\end{pmatrix}
\begin{align}
P(\beta) &= 2(x_{11}\beta_{1}+\cdots+x_{1p}\beta_{p})y_1\\
&+ 2(x_{21}\beta_{1}+\cdots+x_{2p}\beta_{p})y_2\\
&+\cdots\\
&+2 (x_{N1}\beta_{1}+\cdots+x_{Np}\beta_{p})y_N\\
P(\beta) &= \sum_{j=1}^Ny_j\sum_{i=1}^px_{ji}\beta_i
\end{align}
ここで$\beta_1$から$\beta_p$についてそれぞれ微分すると
\frac{\delta P}{\delta \beta_1} = 2(x_{11}y_1+x_{21}y_2+\cdots+x_{N1}y_N) \\
\frac{\delta P}{\delta \beta_2} = 2(x_{12}y_1+x_{22}y_2+\cdots+x_{N2}y_N) \\
\cdots \\
\frac{\delta P}{\delta \beta_p} = 2(x_{1p}y_1+x_{2p}y_2+\cdots+x_{Np}y_N)
ここで、XはN行p列、YはN要素ベクトルであり、$\frac{\delta P}{\delta \beta}$をp要素のベクトルと考えると
\frac{\delta P}{\delta \beta}=2X^Ty\\
次に$Q(\beta)=\beta^TX^TX\beta$を\betaについて微分する
\begin{align}
Q(\beta) &=
\begin{pmatrix}
\beta_{1}\\
\beta_{2}\\
\cdots\\
\beta_{p}\\
\end{pmatrix}^T
\begin{pmatrix}
x_{11} & x_{12} & \cdots & x_{1p}\\
x_{21} & x_{22} & \cdots & x_{2p}\\
\cdots & \cdots & \cdots & \cdots\\
x_{N1} & x_{N2} & \cdots & x_{Np}
\end{pmatrix}^T
\begin{pmatrix}
x_{11} & x_{12} & \cdots & x_{1p}\\
x_{21} & x_{22} & \cdots & x_{2p}\\
\cdots & \cdots & \cdots & \cdots\\
x_{N1} & x_{N2} & \cdots & x_{Np}
\end{pmatrix}
\begin{pmatrix}
\beta_{1}\\
\beta_{2}\\
\cdots\\
\beta_{p}\\
\end{pmatrix}\\
&=
\begin{pmatrix}
\beta_{1} & \cdots & \beta_{p}\\
\end{pmatrix}
\begin{pmatrix}
x_{11} & x_{21} & \cdots & x_{N1}\\
x_{12} & x_{22} & \cdots & x_{N2}\\
\cdots & \cdots & \cdots & \cdots\\
x_{1p} & x_{2p} & \cdots & x_{Np}
\end{pmatrix}
\begin{pmatrix}
x_{11} & x_{12} & \cdots & x_{1p}\\
x_{21} & x_{22} & \cdots & x_{2p}\\
\cdots & \cdots & \cdots & \cdots\\
x_{N1} & x_{N2} & \cdots & x_{Np}
\end{pmatrix}
\begin{pmatrix}
\beta_{1}\\
\beta_{2}\\
\cdots\\
\beta_{p}\\
\end{pmatrix}
\end{align}
$X^TX$について考えると
\begin{align}
X^TX
&=
\begin{pmatrix}
x_{11} & x_{21} & \cdots & x_{N1}\\
x_{12} & x_{22} & \cdots & x_{N2}\\
\cdots & \cdots & \cdots & \cdots\\
x_{1p} & x_{2p} & \cdots & x_{Np}
\end{pmatrix}
\begin{pmatrix}
x_{11} & x_{12} & \cdots & x_{1p}\\
x_{21} & x_{22} & \cdots & x_{2p}\\
\cdots & \cdots & \cdots & \cdots\\
x_{N1} & x_{N2} & \cdots & x_{Np}
\end{pmatrix}\\
&=
\begin{pmatrix}
x_{11}x_{11} + x_{21}x_{21} + \cdots + x_{N1}x_{N1}
& x_{11}x_{12} + x_{21}x_{22} + \cdots + x_{N1}x_{N2}
& \cdots
& x_{11}x_{1p} + x_{21}x_{2p} + \cdots + x_{N1}x_{Np}\\
x_{12}x_{11} + x_{22}x_{21} + \cdots + x_{N2}x_{N1}
& \cdots
& \cdots
& \cdots\\
\cdots
& \cdots
& \cdots
& \cdots\\
x_{1p}x_{11} + x_{2p}x_{21} + \cdots + x_{Np}x_{N1}
& \cdots
& \cdots
& x_{1p}x_{1p} + x_{2p}x_{2p} + \cdots + x_{Np}x_{Np}\\
\end{pmatrix}\\
&=
\begin{pmatrix}
\sum_i^N x_{i1}x_{i1} & \sum_i^N x_{i1}x_{i2} & \cdots & \sum_i^N x_{i1}x_{ip}\\
\sum_i^N x_{i2}x_{i1} & \sum_i^N x_{i2}x_{i2} & \cdots & \sum_i^N x_{i2}x_{ip}\\
\cdots & \cdots & \cdots & \cdots & \\
\sum_i^N x_{ip}x_{i1} & \sum_i^N x_{ip}x_{i2} & \cdots & \sum_i^N x_{ip}x_{ip}\\
\end{pmatrix}
\end{align}
各項は行番号をr、列番号をcとして、$\sum_i^N x_{ir}x_{ic}$と一般化できる。これを$\sum_i^N x_{ir}x_{ic} = X^2_{rc}$と表すとして
X^TX \beta
=
\begin{pmatrix}
X^2_{11}\beta_1 + X^2_{12}\beta_2 + \cdots + X^2_{1p}\beta_p\\
X^2_{21}\beta_1 + X^2_{22}\beta_2 + \cdots + X^2_{2p}\beta_p\\
\cdots\\
X^2_{p1}\beta_1 + X^2_{p2}\beta_2 + \cdots + X^2_{pp}\beta_p
\end{pmatrix}
\begin{align}
Q(\beta)=
\beta^T X^TX \beta
=& (X^2_{11}\beta_1 + X^2_{12}\beta_2 + \cdots + X^2_{1p}\beta_p)\beta_1\\
&+ (X^2_{21}\beta_1 + X^2_{22}\beta_2 + \cdots + X^2_{2p}\beta_p)\beta_2\\
&\cdots\\
&+ (X^2_{p1}\beta_1 + X^2_{p2}\beta_2 + \cdots + X^2_{pp}\beta_p)\beta_p\\
\frac{\delta Q}{\delta \beta_1}
=& 2X^2_{11}\beta_1 + X^2_{12}\beta_2 + X^2_{13}\beta_3 +\cdots + X^2_{1p}\beta_p\\
&+ X^2_{21}\beta_2 + X^2_{31}\beta_3 + \cdots + X^2_{p1}\beta_p
\end{align}
$X^2_{rc}=X^2_{cr}$であることから
\frac{\delta Q}{\delta \beta_1}
= 2(X^2_{11}\beta_1 + X^2_{12}\beta_2 + X^2_{13}\beta_3 +\cdots + X^2_{1p}\beta_p)\\
\frac{\delta Q}{\delta \beta_2}
= 2(X^2_{21}\beta_1 + X^2_{22}\beta_2 + X^2_{23}\beta_3 +\cdots + X^2_{2p}\beta_p)\\
\frac{\delta Q}{\delta \beta_p}
= 2(X^2_{p1}\beta_1 + X^2_{p2}\beta_2 + X^2_{p3}\beta_3 +\cdots + X^2_{pp}\beta_p)\\
これは$X^2\beta=X^TX\beta$である
つまり
\frac{\delta RSS(\beta)}{\delta \beta} = 2X^TX\beta-2X^Ty=0\\
X^TX\beta = X^Ty
となればよく、$X^TX$が正則であるならば、両辺に左からから$(X^TX)^{-1}$をかけて
\beta=(X^TX)^{-1}X^Ty
Rで実装すると
y <- c(80,90,120,110,130,170) #output vector
x0 <- c(1,1,1,1,1,1) #intersept
x1 <- c(5,5,7,5,8,2)
x2 <- c(6,8,10,12,12,12)
X <- cbind(x0, x1, x2) #input matrix
b <- solve(t(X) %*% X) %*% t(X) %*% y
b
lm(y~x1+x2)