はじめに

このノートでは重回帰分析のパラメタを準ニュートン法で計算する方法についてまとめておく。

準ニュートン法

ニュートン法はテイラー展開で2次近似して使用することが多いらしく、2次近似するということは、2階の偏微分が必要になる。偏微分はパラメタが多くなると計算が大変になるため、その点を改善しているのが準ニュートン法(Quasi-Newton Method)とのこと。準ニュートン法はヤコビアンやヘッシアンの近似行列を利用することで、更新を行う。ここではQuasi-Newton method of Broydenの方法をまとめておく。準ニュートン法には他にも多くの近似行列を計算する手法がある。

\[ \begin{eqnarray} s_{t} &=& X_{t+1} - X_{t} \\ y_{t} &=& F_{t+1}(X_{t+1}) - F_{t}(X_{t}) \\ H_{t+1} &=& H_{t} + \left( \frac{s_{t} - H_{t} y_{t}}{s^{T}_{t}H_{t}y_{t}}\right)s^{T}_{t}H_{t} \\ F(X) &=& \begin{pmatrix} f_{1}(x_{1}, \dots, x_{n}) \\ f_{2}(x_{1}, \dots, x_{n}) \\ \vdots \\ f_{n}(x_{1}, \dots, x_{n}) \\ \end{pmatrix} \end{eqnarray} \]

重回帰分析

細かい説明はしないが、重回帰分析は複数の説明変数から目的変数の値を計算する。説明変数を下記の行列で表し、

\[ \begin{eqnarray} \boldsymbol{ X } = \left( \begin{array}{ccc} 1 & \cdots & x_{ 1k } \\ \vdots & \ddots & \vdots \\ 1 & \ldots & x_{ nk } \end{array} \right) \end{eqnarray} \]

目的変数をベクトルとして定める。

\[ \begin{eqnarray} \boldsymbol{ Y } = \left( \begin{array}{c} y_{ 1 } \\ \vdots \\ y_{ n } \end{array} \right) \end{eqnarray} \]

そして、パラメタベクトルを用いて、

\[ \begin{eqnarray} \boldsymbol{ \beta } = \left( \begin{array}{c} \beta_{ 1 } \\ \vdots \\ \beta_{ k } \end{array} \right) \end{eqnarray} \]

下記の線形重回帰式を推定する。\(\boldsymbol{ e }\)は誤差であり、

\[ \boldsymbol{ Y } = \boldsymbol{ X } \boldsymbol{ \beta } + \boldsymbol{ e } \Leftrightarrow \boldsymbol{ e } = \boldsymbol{ Y } - \boldsymbol{ X } \boldsymbol{ \beta } \]

これを最小にするパラメタを求める。この2乗誤差を最小にするために、各パラメタで偏微分する。

\[ \begin{eqnarray} \frac{ \partial \boldsymbol{ e^{\prime} }\boldsymbol{ e } }{ \partial \boldsymbol{ \beta } } &=& \boldsymbol{ 0 } \end{eqnarray} \]

2乗誤差は下記の形に展開でき、

\[ \boldsymbol{ e^{\prime} }\boldsymbol{ e } = \boldsymbol{ y^{\prime} }\boldsymbol{ y } - \boldsymbol{ y^{\prime} }\boldsymbol{ X }\boldsymbol{ \beta } - \boldsymbol{ \beta^{\prime} }\boldsymbol{ X^{\prime} }\boldsymbol{ y } + \boldsymbol{ \beta^{\prime} }\boldsymbol{ X^{\prime} }\boldsymbol{ X }\boldsymbol{ \beta} \]

偏微分すると、

\[ \begin{eqnarray} \frac{ \partial \boldsymbol{ e^{\prime} }\boldsymbol{ e } }{ \partial \boldsymbol{ \beta } } = -2\boldsymbol{ X^{\prime} }\boldsymbol{ y }+ 2\boldsymbol{ X^{\prime} }\boldsymbol{ X }\boldsymbol{ \beta } = \boldsymbol{ 0 } \end{eqnarray} \]

最終的には下記の推定式が得られる。

\[ \begin{eqnarray} -2\boldsymbol{ X^{\prime} }\boldsymbol{ y }+ 2\boldsymbol{ X^{\prime} }\boldsymbol{ X }\boldsymbol{ \beta } = \boldsymbol{ 0 } \Leftrightarrow \boldsymbol{ \beta }=(\boldsymbol{ X^{\prime} }\boldsymbol{ X })^{-1}\boldsymbol{ X^{\prime} }\boldsymbol{ y} \end{eqnarray} \]

重回帰分析と勾配降下法

勾配法で計算したい回帰係数の値を予め確認しておく。

library(palmerpenguins)
library(tidyverse)

penguins_dropna <- penguins %>% drop_na()
Y <- penguins_dropna %>% pull(body_mass_g)
X <- penguins_dropna %>% 
  mutate(x0 = 1) %>% 
  select(x0, x1 = bill_length_mm, x2 = bill_depth_mm) %>%
  as.matrix()

solve(t(X)%*%X)%*%t(X)%*%Y

##          [,1]
## x0 3413.45185
## x1   74.81263
## x2 -145.50718

勾配法では少し時間がかかるが、期待したパラメタが計算できている。

B <- rep(50, ncol(X))
eta <- 10^-6
iter <- 1000000

for(i in 1:iter){
  B <- B - eta * (-2*t(X) %*% Y + 2*t(X) %*% X %*% B)
  if (i %% 100000 == 0) cat(sprintf("%d times: (x0=%2.3f, x1=%2.3f, x2=%2.3f)\n",
                                  i, B[1], B[2], B[3]))
  # 2乗誤差
  # Y_hat <- X %*% beta
  # S <- sum((Y - Y_hat)^2)
}

## 100000 times: (x0=1069.790, x1=99.565, x2=-73.149)
## 200000 times: (x0=1781.445, x1=92.049, x2=-95.121)
## 300000 times: (x0=2277.006, x1=86.815, x2=-110.421)
## 400000 times: (x0=2622.089, x1=83.170, x2=-121.075)
## 500000 times: (x0=2862.387, x1=80.633, x2=-128.494)
## 600000 times: (x0=3029.718, x1=78.865, x2=-133.660)
## 700000 times: (x0=3146.239, x1=77.635, x2=-137.257)
## 800000 times: (x0=3227.379, x1=76.778, x2=-139.762)
## 900000 times: (x0=3283.880, x1=76.181, x2=-141.507)
## 1000000 times: (x0=3323.225, x1=75.766, x2=-142.722)

重回帰分析と準ニュートン法

準ニュートン法は、勾配法よりも収束が早く、期待したパラメタも計算できている。

# Quasi-Newton Algorithm of Broyden
f <- function(b){ c(-2*t(X) %*% Y + 2*t(X) %*% X %*% b) }

iter <- 100000
eta <- 0.001
B <- rep(10, ncol(X))
H <- diag(f(B))

for(i in seq_along(1:iter)){
  # if(sum(abs(f(B))) > 10^(-4)){ # }
    B_pre <- B
    B <- B - eta * H %*% f(B)
    s <- B - B_pre
    y <- f(B) - f(B_pre)
    H <- H + ((s - H %*% y) / as.numeric(t(s) %*% H %*% y)) %*% t(s) %*% H
    if (i %% 10000 == 0) cat(sprintf("%d times: (x0=%2.3f, x1=%2.3f, x2=%2.3f)\n",
                                        i, B[1], B[2], B[3]))
    
    if(sum(abs(f(B))) < 0.0001){
      cat(sprintf("The algorithm converged.\n%d times: (x0=%2.3f, x1=%2.3f, x2=%2.3f)",
                  i, B[1], B[2], B[3]))
      break
      }
}

## 10000 times: (x0=206888927.523, x1=-513884106.841, x2=-74083172.668)
## 20000 times: (x0=12759.163, x1=-23139.055, x2=-3492.085)
## 30000 times: (x0=3413.874, x1=73.764, x2=-145.658)
## 40000 times: (x0=3413.452, x1=74.813, x2=-145.507)
## 50000 times: (x0=3413.452, x1=74.813, x2=-145.507)
## The algorithm converged.
## 53750 times: (x0=3413.452, x1=74.813, x2=-145.507)