重回帰モデルの最適解の解釈：単回帰との整合性

本記事では，「重回帰モデルの最小二乗解は，単回帰モデルと同じ構造をもつ」という主張を示します．
具体的には，デザイン行列 $X$ の列を中心化した行列 $\Delta X$ を導入すると，$\Delta\bold{y}$（平均を引いた後の出力変数）と $\Delta X$（平均を引いた後の入力変数）のあいだに，単回帰で見られる

$$\Delta y = \Delta x\cdot \frac{\mathrm{Cov}(x,y)}{\mathrm{Var}(x)}$$

と同じ形が成立することを示します．

単回帰の最適解の解釈については，先日，こちらの記事（単回帰分析の最適解の解釈）にまとめました．本記事はこの重回帰版で，両者に統一的な理解を与えることを目的としています．

1. 重回帰モデルの基本式と最適解

1.1 重回帰モデルの基本式

重回帰分析では，複数の説明変数を用いて目的変数を予測します．観測データを行列の形で表すと，線形モデルは次のように記述できます．

$$\bold{y} = X\bold{w} + \boldsymbol{\epsilon}$$

ここで，

• $\bold{y} \in \mathbb{R}^{N}$ は目的変数（応答変数）のベクトル
• $X \in \mathbb{R}^{N \times (p+1)}$ はデザイン行列（説明変数を含む行列）
• $\bold{w} \in \mathbb{R}^{(p+1)}$ は回帰係数ベクトル
• $\boldsymbol{\epsilon} \in \mathbb{R}^{N}$ は誤差項（ノイズ）

デザイン行列 $X$ は，通常，$p$ 個の説明変数に加えて，切片（バイアス）を表す $1$ の列を先頭に追加した形をとります．

$$X = \begin{bmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1p} \\ 1 & x_{21} & x_{22} & \cdots & x_{2p} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{N1} & x_{N2} & \cdots & x_{Np} \end{bmatrix}.$$

このモデルでは，$y$ の各要素は次のように表せます．

$$y_i = w_0 + w_1 x_{i1} + w_2 x_{i2} + \cdots + w_p x_{ip} + \epsilon_i.$$

---

1.2 最適解（最小二乗推定量）

(1) 最小二乗法の考え方

回帰係数ベクトル $\bold{w}$ を求めるには，誤差 $\boldsymbol{\epsilon}$ の二乗和（残差平方和）を最小化する最適化問題を解きます．すなわち，目的関数を

$$E(\bold{w}) = \sum_{i=1}^{N} (y_i - X_i \bold{w})^2$$

と定めます．行列形式で表すと，

$$E(\bold{w}) = (\bold{y} - X\bold{w})^T (\bold{y} - X\bold{w}) = \| \bold{y} - X\bold{w} \|_2^2$$

となります．これを最小化するために，$\bold{w}$ に関して偏微分を行います．

(2) 偏微分による最適解の導出

$E(\bold{w})$ を $\bold{w}$ で微分すると，

$$\frac{\partial E(\bold{w})}{\partial \bold{w}} = -2 X^T (\bold{y} - X\bold{w}).$$

これがゼロとなる点が最適解なので，

$$X^T (\bold{y} - X\bold{w}) = \bold{0}$$

を解くことになります．これを整理すると，

$$X^T X \bold{w} = X^T \bold{y}.$$

これが 正規方程式（Normal Equation） です．

(3) 最小二乗推定量

正規方程式を解くと，最適解（最小二乗推定量）は

$$\hat{\bold{w}} = (X^T X)^{-1} X^T y.$$

となります．ただし，$X^T X$ が正則（逆行列を持つ）である場合に限ります．

(4) 予測値の表現

この推定量を用いると，目的変数 $\bold{y}$ の予測値（回帰直線上の値）は

$$\hat{\bold{y}} = X \hat{\bold{w}} = X (X^T X)^{-1} X^T \bold{y}$$

となります．ここで

$$P := X(X^TX)^{-1}X^T$$

を $X$ の射影行列（Projection Matrix）と呼びます．すなわち，$P$ は $\bold{y}$ を $X$ の列空間に射影する行列になっています．

2. 射影行列の説明

2.1 射影行列とは

$X$ の射影行列

$$P = X(X^TX)^{-1}X^T$$

は，$\bold{y}$ を説明変数 $X$ の列空間

$$\mathrm{span}\{ \bold{1},\ \bold{x}_{*1},\ \cdots,\ \bold{x}_{*p} \} = \{X\bold{w}\,|\,\bold{w}\in\mathbb{R}^{p+1}\}$$

に射影する役割を果たします．しかも，$P$ は対称なので， $X$ の列空間への直交射影を与えます．これらのことは，次の二つの事実から理解できます．

(1) $P$ の冪等性

$X$ の射影行列 $P = X(X^TX)^{-1}X^T$ は，ベキ等性

$$P^2 = P$$

を満たします．実際に計算してみると，

$$P^2 = X(X^TX)^{-1}X^T \, X(X^TX)^{-1}X^T = X(X^TX)^{-1}X^T = P$$

となることが確認できます．この事実は，一度射影したものをさらに射影しても変わらないということを意味します．

(2) $P\bold{y}$ が $X$ の列空間への直交射影になっていること

$P\bold{y}$ から $\bold{y}$ に伸びるベクトル（残差ベクトル）

$$\boldsymbol{\epsilon} := \bold{y} - P\bold{y} = (I-P) \bold{y}$$

は $X$ の各列と $\boldsymbol{\epsilon}$ と直交します．すなわち，

$$X^T\boldsymbol{\epsilon} = \begin{bmatrix} \bold{1}^T\boldsymbol{\epsilon} \\ \bold{x}_{*1}^T\boldsymbol{\epsilon} \\ \vdots \\ \bold{x}_{*p}^T\boldsymbol{\epsilon} \end{bmatrix} = \bold{0}$$

が成り立ちます．このことを確認しましょう．

左辺の残差ベクトル $\boldsymbol{\epsilon}$ を展開すると

$$X^T\boldsymbol{\epsilon} = X^T(I-P)\bold{y} = (X^T - X^T P)\bold{y}$$

のような式が現れます．ここで，

$$X^T P = X^TX(X^TX)^{-1}X^T = X^T$$

となることから，

$$X^T \boldsymbol{\epsilon} = \bold{0}$$

が導かれます．こうして，残差ベクトル $\boldsymbol{\epsilon}$ が $X$ の列空間と直交することが分かります．

---

2.2 等価な射影行列

行列 $X$ の列空間は，列ごとに平均を引くなどの列基本変形をしても，張られる空間は変わりません．各説明変数 $\bold{x}_{*j}\ (j=1,\cdots,p)$ の平均

$$\bar{x}_{*j} := \frac{1}{N}\sum_i x_{ij}$$

を引いて中心化した列を

$$\Delta\bold{x}_{*j} := \bold{x}_{*j}-\bar{x}_{*1}\bold{1}$$

と定義し，これを集めた行列を

$$\Delta X = \begin{bmatrix}\begin{array}{c:ccc}\bold{1} & (\bold{x}_{*1}-\bar{x}_{*1}\bold{1}) & \cdots & (\bold{x}_{*p}-\bar{x}_{*p}\bold{1})\end{array}\end{bmatrix} = \begin{bmatrix}\begin{array}{c:ccc}\bold{1} & \Delta\bold{x}_{*1} & \cdots & \Delta\bold{x}_{*p}\end{array}\end{bmatrix}$$

とします．$\Delta X$ の列空間は $X$ の列空間と同じになりますので，射影行列も

$$P = X(X^TX)^{-1}X^T = \Delta X(\Delta X^T\Delta X)^{-1}\Delta X^T$$

という等価な式で書けます．

3. 単回帰と重回帰のアナロジー

3.1 単回帰モデルの復習

説明変数が 1 個の単回帰モデル

$$y = w_o + w_1 x + \epsilon$$

の最適解は

$$\hat{w}_0 = \bar{y} - \hat{w}_1 \bar{x},\quad \hat{w}_1 = \frac{\mathrm{Cov}(x,y)}{\mathrm{Var}(x)} = \frac{\sigma_{xy}}{\sigma_x^2}$$

となります．また，中心化した変数 $\Delta x := x - \bar{x},\ \Delta y := y - \bar{y}$ を用いると，

$$\Delta y = \Delta x \frac{\sigma_{xy}}{\sigma_x^2}$$

という形をしています．

3.2 重回帰モデル

重回帰モデルでは，目的変数 $y$ の予測値は

$$\hat{\bold{y}} = X (X^T X)^{-1} X^T \bold{y}$$

と書けましたが，射影行列の性質から，

$$\hat{\bold{y}} = \Delta X(\Delta X^T\Delta X)^{-1}\Delta X^T \bold{y}$$

と書くことができることを確認しました．この計算を進めてみます．

(1) 行列 $\Delta X^T\Delta X$ の構造

中心化行列 $\Delta X$ は

$$\Delta X = \begin{bmatrix}\begin{array}{c:ccc}\bold{1} & \Delta\bold{x}_{*1} & \cdots & \Delta\bold{x}_{*p}\end{array}\end{bmatrix} = \begin{bmatrix}\begin{array}{c:c}\bold{1} & \Delta\bold{X}_{1:p}\end{array}\end{bmatrix}$$

という形をしており，転置して掛け合わせると

$$\Delta X^T\Delta X = \begin{bmatrix}\begin{array}{c:c} \bold{1}^T\bold{1} & \bold{1}^T \Delta\bold{X}_{1:p} \\ \hdashline (\Delta\bold{X}_{1:p})^T \bold{1} & (\Delta\bold{X}_{1:p})^T \Delta\bold{X}_{1:p} \end{array}\end{bmatrix}.$$

ここで

• $\bold{1}^T \bold{1} = N$.

• $j=1,\cdots,p\ ;\ \bold{1}^T \Delta\bold{x}_{*j} = \sum_i (x_{ij} - \bar{x}_{*j}) = 0$（各列を中心化しているため）．

• $(\Delta\bold{X}_{1:p})^T \Delta\bold{X}_{1:p} = N \sigma_p$，ただし，$\Sigma_p$ は説明変数 $\bold{x}_{1},\cdots\bold{x}_{p}$ の共分散行列

  $$\Sigma_{p} := \begin{bmatrix} \sigma_{11} & \cdots & \sigma_{1p} \\ \vdots & \ddots & \vdots \\ \sigma_{p1} & \cdots & \sigma_{pp} \end{bmatrix}.$$

結局，$\Delta X^T\Delta X$ は

$$\Delta X^T\Delta X = N\begin{bmatrix}\begin{array}{c:c} 1 & \bold{0}^T\\\hdashline \bold{0} & \Sigma_p \end{array}\end{bmatrix}$$

のようなブロック対角行列になり，逆行列は

$$(\Delta X^T\Delta X)^{-1} =\frac{1}{N}\begin{bmatrix}\begin{array}{c:c} 1 & \bold{0}^T\\\hdashline \bold{0} & \Sigma_p^{-1} \end{array}\end{bmatrix}$$

のように計算できます．

(2) ベクトル $\Delta X^T \bold{y}$ の構造

同様に

$$\Delta X^T \bold{y} = \begin{bmatrix}\bold{1}^T\bold{y} \\\hdashline \Delta\bold{X}_{1:p}^T \bold{y} \end{bmatrix} = N \begin{bmatrix} \bar{y} \\\hdashline \boldsymbol{\sigma}_{x_{1:p} y}\end{bmatrix}$$

と書くことができます．ここで

$$\boldsymbol\sigma_{x_{1:p}y} := \begin{bmatrix} \sigma_{x_{1} y} \\ \vdots \\ \sigma_{x_{p} y} \end{bmatrix}$$

は各説明変数と $y$ の共分散を並べたベクトルです．

(3) 予測値 $\hat{\bold{y}}$ の最終式

予測値 $\hat{\bold{y}}$ は射影行列を使って

$$\hat{\bold{y}} = P\bold{y} = \Delta X(\Delta X^T\Delta X)^{-1}\Delta X^T \bold{y}$$

と書けます．ここで，$(\Delta X^T\Delta X)^{-1}$ の係数 $1/N$ と $\Delta X^T y$ の係数 $N$ は互いに打ち消し合って

$$\hat{\bold{y}} = \Delta X\ \begin{bmatrix}\begin{array}{c:c} 1 & \bold{0}^T\\\hdashline \bold{0} & \Sigma_p^{-1} \end{array}\end{bmatrix} \begin{bmatrix} \bar{y} \\\hdashline \bold\sigma_{x_{1:p}y} \end{bmatrix} = \Delta X \begin{bmatrix} \bar{y} \\\hdashline \Sigma_p^{-1} \bold\sigma_{x_{1:p}y} \end{bmatrix} = \bar{\bold{y}}\,\bold{1} + \Delta\bold{X}_{1:p}\, (\Sigma_p^{-1} \bold\sigma_{x_{1:p}y})$$

と書くことができます．また，$\Delta \bold{y} := \bold{y} - \bar{y}\,\bold{1}$ を用いると，予測値 $\hat{\bold{y}}$ は

$$\Delta \bold{y} = \Delta\bold{X}_{1:p}\, (\Sigma_p^{-1} \bold\sigma_{x_{1:p}y})$$

を満たすことになります．ここから次の2つのことが分かります．

• 入力データの平均ベクトル $(\bar{x}_{*1},\cdots,\bar{x}_{*1})$ に出力の平均 $\bar{y}$ が対応する
• デザイン行列の 1 列目 $\bold{1}$ はデータを中心化することに寄与する

4. まとめ

単回帰と重回帰の最適解を並べると

• 単回帰：$\Delta y = \Delta x \frac{\sigma_{xy}}{\sigma_x^2}$
• 重回帰：$\Delta \bold{y} = \Delta\bold{X}_{1:p}\, (\Sigma_p^{-1} \boldsymbol{\sigma}_{x_{1:p}y})$

のように，同じ形をしていることが分かります．どちらも「データを平均からの差分に変換してみると，変数の共分散や分散の比によって予測値が決まる」という点で統一的に理解できます．