線形回帰と最小二乗法：正規方程式を直交射影として読む

0. この記事の要点

測定した点が何十個もあると、それらをすべて通る直線は引けません。そこで「各点とのずれを二乗して合計した値が一番小さい直線」を選ぶのが最小二乗法です。この記事では、その直線が連立一次方程式を解くだけで求まること、そして計算の意味が「平面に垂線を下ろす」という図で見えてくることを説明します。以下の要点は内容を数学の言葉で並べたものなので、分からない言葉があっても、後の節で順に説明します。

• 線形回帰の学習とは、観測ベクトル $\boldsymbol{y}$ に対して $\boldsymbol{y} \approx X\boldsymbol{w}$ をできるだけよく満たす $\boldsymbol{w}$ を探すことです。切片 $b$ は、計画行列 $X$ に全成分 $1$ の列を追加することで重みの一部として吸収できます。
• 二乗和誤差 $L(\boldsymbol{w}) = \lVert \boldsymbol{y} - X\boldsymbol{w}\rVert^2$ の最小化は、微分を経由しなくても、正規方程式 $X^{\mathsf{T}}X\boldsymbol{w} = X^{\mathsf{T}}\boldsymbol{y}$ を解くことと完全に同値です（定理 3.3）。学習が連立一次方程式に化けるところが、このモデルの要点です。
• 正規方程式の解 $\hat{\boldsymbol{w}}$ が与える予測 $X\hat{\boldsymbol{w}}$ は、$\boldsymbol{y}$ の列空間 $\operatorname{Im} X$（$X$ の列の一次結合全体）への直交射影（垂線を下ろして得られる点）にほかなりません（定理 4.2）。残差はすべての説明変数と直交します。
• 正規方程式は $X$ や $\boldsymbol{y}$ が何であっても必ず解を持ちます。解が一意になるのは $\operatorname{rank} X = p$（列がすべて一次独立）のときで、そのとき $\hat{\boldsymbol{w}} = (X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}\boldsymbol{y}$ です。ランクが落ちても予測 $X\hat{\boldsymbol{w}}$ のほうは一意に決まります。
• 実装では $(X^{\mathsf{T}}X)^{-1}$ を陽に計算しないでください。$\kappa_2(X^{\mathsf{T}}X) = \kappa_2(X)^2$ と条件数（計算の誤差の増えやすさの目安）が二乗されるため、QR 分解や特異値分解（どちらも行列を扱いやすい行列の積に分ける方法）を経由するほうが安全です。

この記事で使う事前知識

• 部分空間: ベクトルの和とスカラー倍で閉じている集合のこと。列空間 $\operatorname{Im}X$ や核 $\operatorname{Ker}X$ がその例です。→ 部分空間
• 一次独立: どのベクトルも、残りのベクトルの一次結合では表せない、という性質。→ 線形独立・基底
• 次元: 空間の中で一次独立に取れるベクトルの最大個数。ランクはこの次元の考え方と結びつきます。→ 次元
• 凸関数: グラフの 2 点を結ぶ線分が、グラフより上にある関数。停留点が最小点になります。→ 凸関数
• 経験リスク最小化: 訓練データでの損失の平均を最小にしてパラメータを決める、という学習の見方。→ なぜ機械学習に数学が必要か

1. 動機

1.1. 情報が多すぎるという問題

中学校で習ったとおり、平面上の相異なる 2 点を通る直線はただ 1 本に決まります。では、測定誤差を含む 100 点が与えられたとき、それらを「通る」直線はどれでしょうか。答えは「1 本もない」です。未知数は傾きと切片の 2 個しかないのに、方程式は 100 本あります。こうした過剰決定系（方程式の数が未知数の数より多い連立一次方程式）は、普通の意味では解を持ちません。

この状況は、19 世紀初頭の天文学と測地学ではごく日常的なものでした。惑星や彗星の軌道要素は 6 個ですが、観測は何十回も行われます。当時の実務家は、観測をいくつか選んで方程式を立てたり、観測を組にして平均したりと、場当たり的な処理をしていました。「すべての観測を公平に使い、しかも一意な答えを返す手続き」が求められていたのです。

その手続きが最小二乗法です。1801 年 1 月 1 日にジュゼッペ・ピアッツィが小惑星ケレスを発見しましたが、41 日間の観測ののちに太陽の方向へ入って見失われました。当時 24 歳のガウスは、わずかな観測データから軌道を計算して再発見の位置を予告し、同じ年の暮れに天文学者たちがその予告どおりの場所でケレスを再発見します。この事件は最小二乗法の威力を広く知らしめました。方法自体は 1805 年にルジャンドルが著書『Nouvelles méthodes pour la détermination des orbites des comètes』の付録で “méthode des moindres carrés” として公表し、ガウスは 1809 年の『Theoria Motus』で「自分は 1795 年から使っていた」と主張して、有名な優先権論争になりました。

1.2. なぜ「二乗」なのか

誤差の測り方は二乗和だけではありません。絶対値の和 $\sum_i |y_i - \hat{y}_i|$ でも、最大誤差 $\max_i |y_i - \hat{y}_i|$ でもよいはずです。それでも二乗和が標準になったのには、はっきりした理由があります。

1. 微分できます。 絶対値は原点で微分できませんが、二乗は至るところ滑らかで、しかも $\boldsymbol{w}$ の二次式です。二次式の停留点を求める条件は一次方程式になります。
2. 幾何と直結します。 二乗和はユークリッドノルムの二乗、つまり内積から来る距離です。したがって「誤差を最小にする」は「部分空間へ垂線を下ろす」と同じ意味になり、ピタゴラスの定理がそのまま使えます（定理 4.2）。
3. 確率的な意味があります。 観測ノイズが独立同分布の正規分布に従うと仮定すると、対数尤度の最大化が二乗和誤差の最小化とぴったり一致します。これはガウス自身の正当化であり、詳しくは 確率論とベイズ統計の役割 の ガウス雑音の下で最尤推定は最小二乗法（命題 3.3）［確率論とベイズ統計］ で扱います。

一方で弱点もあります。二乗は大きな誤差をさらに強調するので、外れ値 1 点が解を大きく動かします。この弱点への対処（ロバスト回帰、正則化）は本記事の後半と演習で触れます。

1.3. 機械学習における位置づけ

なぜ機械学習に数学が必要か？ で見たとおり、機械学習のアルゴリズムは数学の言葉で書かれると急に見通しがよくなります。線形回帰は、教師あり学習のなかで学習問題が閉じた式で解ける唯一といってよいモデルです。ロジスティック回帰やニューラルネットワークでは、損失関数が $\boldsymbol{w}$ の二次式でなくなるため閉じた解が消え、勾配降下法 のような反復解法に頼ることになります。だからこそ、閉じた解が存在するこの場合に「何が起きているのか」を完全に理解しておく価値があります。ここで現れる直交射影・ランク条件・条件数という三つ組は、主成分分析 にも ロジスティック回帰 の反復重み付き最小二乗にも、そのまま持ち越されます。

2. 準備：データを行列に並べる

データは $n$ 組の観測 $(\boldsymbol{x}_i, y_i)$（$i = 1, \dots, n$）とします。$\boldsymbol{x}_i \in \mathbb{R}^d$ は説明変数（特徴量）のベクトル、$y_i \in \mathbb{R}$ は目的変数です。予測したい関係は

$$\hat{y} = w_1 x_1 + w_2 x_2 + \cdots + w_d x_d + b$$

という形をしています。$b$ は切片（バイアス項）です。この $d+1$ 個の未知数をまとめて扱うために、次の記法を用意します。

定義 2.1（計画行列と線形回帰モデル）↗

$n$ 個の観測 $(\boldsymbol{x}_i, y_i)$、$\boldsymbol{x}_i = (x_{i1}, \dots, x_{id})^{\mathsf{T}} \in \mathbb{R}^d$、$y_i \in \mathbb{R}$ に対し、$p = d + 1$ とおき

$$X = \begin{pmatrix} 1 & x_{11} & \cdots & x_{1d} \\ 1 & x_{21} & \cdots & x_{2d} \\ \vdots & \vdots & & \vdots \\ 1 & x_{n1} & \cdots & x_{nd} \end{pmatrix} \in \mathbb{R}^{n \times p}, \qquad \boldsymbol{y} = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix} \in \mathbb{R}^{n}$$

とする。$X$ を計画行列（design matrix）という。第 1 列は全成分が $1$ のベクトル $\boldsymbol{1}$ である。パラメータを $\boldsymbol{w} = (b, w_1, \dots, w_d)^{\mathsf{T}} \in \mathbb{R}^{p}$ とすると、$n$ 個の予測値はまとめて $X\boldsymbol{w} \in \mathbb{R}^n$ と書ける。これを線形回帰モデルという。

第 1 列に $\boldsymbol{1}$ を置いたことで、切片は「定数 $1$ という特徴量に対する重み」になりました。以後、切片を特別扱いする必要はありません。$X$ の第 $i$ 行は $i$ 番目のデータ、第 $j$ 列は $j$ 番目の特徴量が全データにわたって並んだベクトルです。行はデータ、列は特徴量という対応を、以下ずっと使います。

ここで「線形」の意味を正確にしておきます。線形回帰の「線形」は、説明変数についてではなくパラメータ $\boldsymbol{w}$ について線形という意味です。この区別が、線形回帰の応用範囲を大きく広げます。

定義 2.2（特徴写像による拡張）↗

写像 $\varphi \colon \mathbb{R}^d \to \mathbb{R}^{p}$ を任意に固定し、計画行列の第 $i$ 行を $\varphi(\boldsymbol{x}_i)^{\mathsf{T}}$ で置き換えたものを $\Phi \in \mathbb{R}^{n \times p}$ とする。モデル $\hat{y} = \boldsymbol{w}^{\mathsf{T}}\varphi(\boldsymbol{x})$ を、特徴写像 $\varphi$ による線形モデルという。$\varphi$ がどれほど非線形でも、$\boldsymbol{w}$ に関しては線形である。

たとえば $d = 1$ で $\varphi(x) = (1, x, x^2, x^3)^{\mathsf{T}}$ とすれば 3 次多項式の当てはめになりますが、未知数 $\boldsymbol{w}$ については依然として線形です。したがって以下で作る道具立ては、そのまま多項式回帰・三角関数展開・動径基底関数モデルに適用できます。実際の数値例は 例 6.2 で見ます。

ノート

以下、$\mathbb{R}^n$ には標準内積 $\langle \boldsymbol{u}, \boldsymbol{v}\rangle = \boldsymbol{u}^{\mathsf{T}}\boldsymbol{v} = \sum_{i} u_i v_i$ とノルム $\lVert \boldsymbol{u}\rVert = \sqrt{\langle \boldsymbol{u}, \boldsymbol{u}\rangle}$ を入れます。行列 $A$ に対し $\operatorname{Im} A = \{A\boldsymbol{v} : \boldsymbol{v}\}$（列空間、値域）、$\operatorname{Ker} A = \{\boldsymbol{v} : A\boldsymbol{v} = \boldsymbol{0}\}$（核）と書きます。線形代数からは、次元定理 $\operatorname{rank} A + \dim \operatorname{Ker} A = (\text{列数})$、行階数と列階数の一致 $\operatorname{rank} A = \operatorname{rank} A^{\mathsf{T}}$、および部分空間 $W \subseteq \mathbb{R}^n$ に対する直交分解 $\mathbb{R}^n = W \oplus W^{\perp}$ を既知として使います。前二者は 行列と連立一次方程式 の 次元定理（定理 8.3）［行列と連立一次方程式］ と 行階数と列階数の一致（定理 8.2）［行列と連立一次方程式］、最後のものは 内積空間とグラム・シュミット直交化 の 直交分解定理（定理 7.1）［内積空間とグラム・シュミット直交化］ を参照してください。

3. 最小二乗問題と正規方程式

3.1. 問題の定式化

$n > p$ のとき、$X\boldsymbol{w} = \boldsymbol{y}$ をぴったり満たす $\boldsymbol{w}$ は一般に存在しません。そこで「ぴったり」をあきらめ、ずれの大きさを最小にする $\boldsymbol{w}$ を選びます。

定義 3.1（最小二乗問題）↗

$X \in \mathbb{R}^{n \times p}$、$\boldsymbol{y} \in \mathbb{R}^n$ に対し、残差ベクトルを $\boldsymbol{r}(\boldsymbol{w}) = \boldsymbol{y} - X\boldsymbol{w}$、二乗和誤差（残差平方和）を

$$L(\boldsymbol{w}) = \lVert \boldsymbol{y} - X\boldsymbol{w}\rVert^2 = \sum_{i=1}^{n} \bigl(y_i - (X\boldsymbol{w})_i\bigr)^2$$

と定める。$L$ を $\mathbb{R}^p$ 全体で最小にする $\boldsymbol{w}$、すなわち

$$\hat{\boldsymbol{w}} \in \operatorname*{arg\,min}_{\boldsymbol{w} \in \mathbb{R}^{p}} \lVert \boldsymbol{y} - X\boldsymbol{w}\rVert^2$$

を満たす $\hat{\boldsymbol{w}}$ を最小二乗解という。

$L$ は $\boldsymbol{w}$ の各成分についての二次式であり、下に有界（つねに $L \ge 0$）です。以下で見るように、最小値は必ず達成されます。しかし最小点が 1 つとは限りません。この「存在するが一意とは限らない」という構造を正確に捉えることが、この節と次節の目標です。

3.2. 二次形式の勾配

まず、あとで使う微分公式を独立に証明しておきます。ここでの偏微分は多変数の意味ですから、必要なら 多変数関数の微分と偏微分 の 偏微分の定義（定義 3.1）［多変数関数の微分と偏微分］ を復習してください。

補題 3.2（二次形式の勾配）✓

$A \in \mathbb{R}^{p \times p}$ を対称行列（$A^{\mathsf{T}} = A$）、$\boldsymbol{b} \in \mathbb{R}^p$、$c \in \mathbb{R}$ とし、関数 $f \colon \mathbb{R}^p \to \mathbb{R}$ を

$$f(\boldsymbol{w}) = \boldsymbol{w}^{\mathsf{T}} A \boldsymbol{w} - 2\boldsymbol{b}^{\mathsf{T}}\boldsymbol{w} + c$$

で定める。このとき $f$ は $C^{\infty}$ 級であり、勾配とヘッセ行列は

$$\nabla f(\boldsymbol{w}) = 2A\boldsymbol{w} - 2\boldsymbol{b}, \qquad \nabla^2 f(\boldsymbol{w}) = 2A$$

で与えられる。

証明（補題 3.2）

成分で書くと $f(\boldsymbol{w}) = \sum_{i=1}^{p}\sum_{j=1}^{p} a_{ij} w_i w_j - 2\sum_{i=1}^{p} b_i w_i + c$ であり、これは $w_1, \dots, w_p$ の多項式ですから $C^{\infty}$ 級です。第 $k$ 成分で偏微分します。積 $w_i w_j$ が $w_k$ を含むのは $i = k$ の場合と $j = k$ の場合で、$i = j = k$ の項 $a_{kk}w_k^2$ は微分すると $2a_{kk}w_k$ となり、両方の数え方の和と一致します。したがって

$$\frac{\partial f}{\partial w_k}(\boldsymbol{w}) = \sum_{j=1}^{p} a_{kj} w_j + \sum_{i=1}^{p} a_{ik} w_i - 2 b_k = (A\boldsymbol{w})_k + (A^{\mathsf{T}}\boldsymbol{w})_k - 2b_k$$

です。仮定 $A^{\mathsf{T}} = A$ より右辺は $2(A\boldsymbol{w})_k - 2b_k$ となり、第一式が出ます。さらにこれを $w_l$ で偏微分すると $\partial^2 f / \partial w_l \partial w_k = 2a_{kl}$ であり、ヘッセ行列は $2A$ です。

∎

3.3. 正規方程式

いよいよ中心となる定理です。微分を一切使わずに証明できることに注意してください。証明が代数的な恒等式だけで済むので、最小値であること（停留点にすぎないのではないこと）が同時に示せます。

定理 3.3（正規方程式）✓

$X \in \mathbb{R}^{n\times p}$、$\boldsymbol{y} \in \mathbb{R}^n$ とし、$L(\boldsymbol{w}) = \lVert \boldsymbol{y} - X\boldsymbol{w}\rVert^2$ とする。$\hat{\boldsymbol{w}} \in \mathbb{R}^p$ について、次の 2 条件は同値である。

1. $\hat{\boldsymbol{w}}$ は $L$ を $\mathbb{R}^p$ 上で最小にする、すなわちすべての $\boldsymbol{w} \in \mathbb{R}^p$ に対し $L(\hat{\boldsymbol{w}}) \le L(\boldsymbol{w})$。
2. $\hat{\boldsymbol{w}}$ は正規方程式 $X^{\mathsf{T}}X\hat{\boldsymbol{w}} = X^{\mathsf{T}}\boldsymbol{y}$ を満たす。

さらにこのとき、任意の $\boldsymbol{w} \in \mathbb{R}^p$ に対して恒等式

$$L(\boldsymbol{w}) = L(\hat{\boldsymbol{w}}) + \lVert X(\boldsymbol{w} - \hat{\boldsymbol{w}})\rVert^2$$

が成り立ち、最小二乗解の全体は $\hat{\boldsymbol{w}} + \operatorname{Ker} X = \{\hat{\boldsymbol{w}} + \boldsymbol{v} : X\boldsymbol{v} = \boldsymbol{0}\}$ と一致する。

証明の方針は次のとおりです。最小点の候補 $\hat{\boldsymbol{w}}$ から $\boldsymbol{v}$ だけ動かしたときの誤差の変化を、$\boldsymbol{v}$ の一次の項と二次の項に分けて書き出します。一次の項の係数が $\boldsymbol{0}$ なら二次の項（$0$ 以上）だけが残るので最小点です。$\boldsymbol{0}$ でなければ、その係数の向きに十分小さく動くと一次の項が二次の項に勝って誤差が減る、というのが後半で $t\boldsymbol{g}$ という動かし方を選ぶ理由です。

証明（定理 3.3）

任意の $\hat{\boldsymbol{w}}, \boldsymbol{v} \in \mathbb{R}^p$ に対し、$\hat{\boldsymbol{r}} = \boldsymbol{y} - X\hat{\boldsymbol{w}}$ とおいて $L(\hat{\boldsymbol{w}} + \boldsymbol{v})$ を展開します。ノルムの二乗は内積なので

$$\begin{aligned} L(\hat{\boldsymbol{w}} + \boldsymbol{v}) &= \lVert \hat{\boldsymbol{r}} - X\boldsymbol{v}\rVert^2 \\ &= \langle \hat{\boldsymbol{r}} - X\boldsymbol{v},\ \hat{\boldsymbol{r}} - X\boldsymbol{v}\rangle \\ &= \lVert \hat{\boldsymbol{r}}\rVert^2 - 2\langle \hat{\boldsymbol{r}},\ X\boldsymbol{v}\rangle + \lVert X\boldsymbol{v}\rVert^2 \end{aligned}$$

となります。ここで中央の項は転置の性質 $\langle \boldsymbol{u}, A\boldsymbol{v}\rangle = \langle A^{\mathsf{T}}\boldsymbol{u}, \boldsymbol{v}\rangle$ により

$$\langle \hat{\boldsymbol{r}},\ X\boldsymbol{v}\rangle = \langle X^{\mathsf{T}}\hat{\boldsymbol{r}},\ \boldsymbol{v}\rangle = \boldsymbol{v}^{\mathsf{T}}\bigl(X^{\mathsf{T}}\boldsymbol{y} - X^{\mathsf{T}}X\hat{\boldsymbol{w}}\bigr)$$

と書けます。そこで $\boldsymbol{g} = X^{\mathsf{T}}\boldsymbol{y} - X^{\mathsf{T}}X\hat{\boldsymbol{w}} = X^{\mathsf{T}}\hat{\boldsymbol{r}}$ とおくと、すべての $\boldsymbol{v}$ について

$$L(\hat{\boldsymbol{w}} + \boldsymbol{v}) = L(\hat{\boldsymbol{w}}) - 2\boldsymbol{v}^{\mathsf{T}}\boldsymbol{g} + \lVert X\boldsymbol{v}\rVert^2$$

が成り立ちます。この 1 本の恒等式から両方向が出ます。

(2) $\Rightarrow$ (1)。 正規方程式が成り立てば $\boldsymbol{g} = \boldsymbol{0}$ ですから、上の恒等式は $L(\hat{\boldsymbol{w}} + \boldsymbol{v}) = L(\hat{\boldsymbol{w}}) + \lVert X\boldsymbol{v}\rVert^2$ になります。$\lVert X\boldsymbol{v}\rVert^2 \ge 0$ なので、任意の $\boldsymbol{v}$ に対し $L(\hat{\boldsymbol{w}} + \boldsymbol{v}) \ge L(\hat{\boldsymbol{w}})$、すなわち $\hat{\boldsymbol{w}}$ は最小点です。$\boldsymbol{w} = \hat{\boldsymbol{w}} + \boldsymbol{v}$ と置き直せば主張の恒等式そのものになります。また等号成立は $\lVert X\boldsymbol{v}\rVert = 0$、つまり $X\boldsymbol{v} = \boldsymbol{0}$ と同値なので、最小点の全体は $\hat{\boldsymbol{w}} + \operatorname{Ker} X$ です。

(1) $\Rightarrow$ (2)。 対偶を示します。$\boldsymbol{g} \ne \boldsymbol{0}$ と仮定し、$\boldsymbol{v} = t\boldsymbol{g}$（$t > 0$）を代入すると

$$L(\hat{\boldsymbol{w}} + t\boldsymbol{g}) - L(\hat{\boldsymbol{w}}) = -2t\lVert \boldsymbol{g}\rVert^2 + t^2 \lVert X\boldsymbol{g}\rVert^2 = t\bigl(t\lVert X\boldsymbol{g}\rVert^2 - 2\lVert \boldsymbol{g}\rVert^2\bigr)$$

です。$X\boldsymbol{g} = \boldsymbol{0}$ ならば右辺は任意の $t > 0$ に対し $-2t\lVert\boldsymbol{g}\rVert^2 < 0$ です。$X\boldsymbol{g} \ne \boldsymbol{0}$ ならば、たとえば $t = \lVert \boldsymbol{g}\rVert^2 / \lVert X\boldsymbol{g}\rVert^2 > 0$ と取れば右辺は $-\lVert\boldsymbol{g}\rVert^4/\lVert X\boldsymbol{g}\rVert^2 < 0$ です。いずれの場合も $L(\hat{\boldsymbol{w}} + t\boldsymbol{g}) < L(\hat{\boldsymbol{w}})$ となり、$\hat{\boldsymbol{w}}$ は最小点ではありません。したがって最小点なら $\boldsymbol{g} = \boldsymbol{0}$、すなわち正規方程式が成り立ちます。

∎

正規方程式 $X^{\mathsf{T}}X\boldsymbol{w} = X^{\mathsf{T}}\boldsymbol{y}$ は $p$ 元の連立一次方程式です。$n$ がどれほど大きくても（データが何億件でも）、解くべき方程式の大きさは特徴量の個数 $p$ だけで決まります。「学習」という言葉のうしろで実際に起きているのは、$p \times p$ の連立一次方程式を解くことなのです。

注意 3.4（偏微分による導出）

微分を使う導き方も確認しておきます。展開すると

$$L(\boldsymbol{w}) = \lVert \boldsymbol{y}\rVert^2 - 2\boldsymbol{y}^{\mathsf{T}}X\boldsymbol{w} + \boldsymbol{w}^{\mathsf{T}}X^{\mathsf{T}}X\boldsymbol{w}$$

であり、$A = X^{\mathsf{T}}X$（$(X^{\mathsf{T}}X)^{\mathsf{T}} = X^{\mathsf{T}}X$ だから対称）、$\boldsymbol{b} = X^{\mathsf{T}}\boldsymbol{y}$、$c = \lVert\boldsymbol{y}\rVert^2$ として 補題 3.2 を適用すれば

$$\nabla L(\boldsymbol{w}) = 2X^{\mathsf{T}}X\boldsymbol{w} - 2X^{\mathsf{T}}\boldsymbol{y}$$

を得ます。$\nabla L(\hat{\boldsymbol{w}}) = \boldsymbol{0}$ が正規方程式です。ただし停留点であることは一般には最小であることを意味しません。ここではヘッセ行列が $2X^{\mathsf{T}}X$ で、任意の $\boldsymbol{v}$ に対し $\boldsymbol{v}^{\mathsf{T}}(X^{\mathsf{T}}X)\boldsymbol{v} = \lVert X\boldsymbol{v}\rVert^2 \ge 0$、すなわち半正定値なので $L$ は凸関数であり、停留点が大域的最小点になります。定理 3.3 の証明は、この凸性の議論を恒等式 1 本に置き換えたものだと見ることができます。

4. 幾何：最小二乗解は直交射影である

4.1. 残差はすべての特徴量と直交する

正規方程式を $X^{\mathsf{T}}(\boldsymbol{y} - X\hat{\boldsymbol{w}}) = \boldsymbol{0}$、つまり $X^{\mathsf{T}}\hat{\boldsymbol{r}} = \boldsymbol{0}$ と書き直してみます。$X^{\mathsf{T}}\hat{\boldsymbol{r}}$ の第 $j$ 成分は $X$ の第 $j$ 列と $\hat{\boldsymbol{r}}$ の内積ですから、この式は

残差ベクトルは、$X$ のすべての列（＝すべての特徴量）と直交する

と読めます。統計的に言えば「残差の中には、手持ちの特徴量の一次結合で説明できる成分がもう残っていない」ということです。最小二乗法は、説明変数で説明しつくせる分をすべて使い切ったところで止まる手続きなのです。この読み方を定理の形にします。

まず 2 次元の例で図を見ておきます。原点を通る直線（1 本の列がつくる「空間」）と、その外にある点 $\boldsymbol{y}$ があるとき、$\boldsymbol{y}$ に一番近い直線上の点は、$\boldsymbol{y}$ から垂線を下ろした足です。このとき垂線（残差）は直線と直交します。列が 2 本なら直線が平面に変わるだけで、同じことが成り立ちます。

2 次元での直交射影。直線上の点のうち y に最も近いのは垂線の足 ŷ で、残差 r = y − ŷ は直線と直交する。

次の補題は、「残差が列と直交する」ことを「$X^{\mathsf{T}}$ をかけると $\boldsymbol{0}$」という式に言い換えるものです。証明では、「すべての $\boldsymbol{v}$ で内積が $0$」という強い条件を使い、$\boldsymbol{v}$ に最も都合のよいベクトル $X^{\mathsf{T}}\boldsymbol{u}$ 自身を入れます。すると内積が長さの二乗になり、それが $0$ なら零ベクトルだと分かります。

補題 4.1（列空間の直交補空間）✓

任意の $X \in \mathbb{R}^{n\times p}$ に対し $(\operatorname{Im} X)^{\perp} = \operatorname{Ker} X^{\mathsf{T}}$ が成り立つ。

証明（補題 4.1）

$\boldsymbol{u} \in (\operatorname{Im}X)^{\perp}$ であることは、すべての $\boldsymbol{v}\in\mathbb{R}^p$ に対して $\langle \boldsymbol{u}, X\boldsymbol{v}\rangle = 0$ であることです。転置の性質より $\langle \boldsymbol{u}, X\boldsymbol{v}\rangle = \langle X^{\mathsf{T}}\boldsymbol{u}, \boldsymbol{v}\rangle$ ですから、これは「すべての $\boldsymbol{v}$ に対し $\langle X^{\mathsf{T}}\boldsymbol{u}, \boldsymbol{v}\rangle = 0$」と同値です。とくに $\boldsymbol{v} = X^{\mathsf{T}}\boldsymbol{u}$ と取れば $\lVert X^{\mathsf{T}}\boldsymbol{u}\rVert^2 = 0$ すなわち $X^{\mathsf{T}}\boldsymbol{u} = \boldsymbol{0}$ を得ます。逆に $X^{\mathsf{T}}\boldsymbol{u} = \boldsymbol{0}$ ならすべての $\boldsymbol{v}$ で内積は $0$ です。よって両者は同値で、集合として一致します。

∎

定理 4.2（最小二乗解の存在と射影表示）✓

$X \in \mathbb{R}^{n\times p}$、$\boldsymbol{y}\in\mathbb{R}^n$ を任意とする。このとき次が成り立つ。

1. 正規方程式 $X^{\mathsf{T}}X\boldsymbol{w} = X^{\mathsf{T}}\boldsymbol{y}$ は少なくとも 1 つの解を持つ。とくに最小二乗解は必ず存在する。
2. $\hat{\boldsymbol{w}}$ を任意の最小二乗解とするとき、$X\hat{\boldsymbol{w}}$ は $\boldsymbol{y}$ の部分空間 $\operatorname{Im}X$ への直交射影に一致する。とくに $X\hat{\boldsymbol{w}}$ は最小二乗解の取り方によらず一意に定まる。
3. $\hat{\boldsymbol{y}} = X\hat{\boldsymbol{w}}$ は、$\operatorname{Im}X$ の元のうち $\boldsymbol{y}$ に最も近い唯一の点である。すなわち $\boldsymbol{z} \in \operatorname{Im}X$、$\boldsymbol{z} \ne \hat{\boldsymbol{y}}$ ならば $\lVert \boldsymbol{y}-\boldsymbol{z}\rVert > \lVert \boldsymbol{y}-\hat{\boldsymbol{y}}\rVert$。

証明（定理 4.2）

(1)。 $W = \operatorname{Im}X$ は $\mathbb{R}^n$ の部分空間ですから、直交分解 $\mathbb{R}^n = W \oplus W^{\perp}$ により $\boldsymbol{y} = \hat{\boldsymbol{y}} + \boldsymbol{s}$（$\hat{\boldsymbol{y}} \in W$、$\boldsymbol{s}\in W^{\perp}$）と一意に分解できます。$\hat{\boldsymbol{y}} \in \operatorname{Im}X$ なので $\hat{\boldsymbol{y}} = X\hat{\boldsymbol{w}}$ となる $\hat{\boldsymbol{w}}$ が存在します。補題 4.1 より $\boldsymbol{s} \in W^{\perp} = \operatorname{Ker}X^{\mathsf{T}}$、すなわち $X^{\mathsf{T}}\boldsymbol{s} = \boldsymbol{0}$ です。したがって

$$X^{\mathsf{T}}X\hat{\boldsymbol{w}} = X^{\mathsf{T}}\hat{\boldsymbol{y}} = X^{\mathsf{T}}(\boldsymbol{y} - \boldsymbol{s}) = X^{\mathsf{T}}\boldsymbol{y}$$

となり、この $\hat{\boldsymbol{w}}$ は正規方程式の解です。定理 3.3 よりこれは最小二乗解です。

(2)。 $\hat{\boldsymbol{w}}$ を任意の最小二乗解とすると、定理 3.3 より $X^{\mathsf{T}}(\boldsymbol{y} - X\hat{\boldsymbol{w}}) = \boldsymbol{0}$、すなわち 補題 4.1 により $\boldsymbol{y} - X\hat{\boldsymbol{w}} \in W^{\perp}$ です。一方 $X\hat{\boldsymbol{w}} \in W$ ですから、$\boldsymbol{y} = X\hat{\boldsymbol{w}} + (\boldsymbol{y}-X\hat{\boldsymbol{w}})$ は $W \oplus W^{\perp}$ に沿った分解になっています。この分解は一意なので、$X\hat{\boldsymbol{w}}$ は (1) で作った $\hat{\boldsymbol{y}}$（＝ $\boldsymbol{y}$ の $W$ への直交射影）に等しく、最小二乗解の選び方に依存しません。

(3)。 $\boldsymbol{z}\in W$ とすると $\hat{\boldsymbol{y}} - \boldsymbol{z} \in W$ であり、$\boldsymbol{y}-\hat{\boldsymbol{y}} \in W^{\perp}$ なので両者は直交します。よって ピタゴラスの定理（注意 4.5）［内積空間とグラム・シュミット直交化］ より

$$\lVert \boldsymbol{y}-\boldsymbol{z}\rVert^2 = \lVert (\boldsymbol{y}-\hat{\boldsymbol{y}}) + (\hat{\boldsymbol{y}}-\boldsymbol{z})\rVert^2 = \lVert \boldsymbol{y}-\hat{\boldsymbol{y}}\rVert^2 + \lVert \hat{\boldsymbol{y}}-\boldsymbol{z}\rVert^2$$

です。$\boldsymbol{z}\ne\hat{\boldsymbol{y}}$ なら第 2 項は正なので、狭義の不等式が成り立ちます。

∎

最小二乗解の幾何。観測ベクトル y は一般に列空間 Im X の外にあり、予測ベクトルは y から Im X へ下ろした垂線の足（直交射影）である。残差 r は列空間全体と直交する。

4.2. ハット行列と決定係数

$\operatorname{rank}X = p$ の場合には、射影を行列の形で書き下せます。

命題 4.3（ハット行列）✓

$X \in \mathbb{R}^{n\times p}$ が $\operatorname{rank}X = p$ を満たすとする（このとき 系 5.2 により $X^{\mathsf{T}}X$ は正則）。

$$P = X(X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}} \in \mathbb{R}^{n\times n}$$

とおくと、次が成り立つ。

1. $P^{\mathsf{T}} = P$ かつ $P^2 = P$（対称冪等）。
2. 任意の $\boldsymbol{y}\in\mathbb{R}^n$ に対し $P\boldsymbol{y}$ は $\boldsymbol{y}$ の $\operatorname{Im}X$ への直交射影であり、$P\boldsymbol{y} = X\hat{\boldsymbol{w}}$。
3. $I_n - P$ は $(\operatorname{Im}X)^{\perp} = \operatorname{Ker}X^{\mathsf{T}}$ への直交射影であり、残差は $\hat{\boldsymbol{r}} = (I_n - P)\boldsymbol{y}$。
4. $\operatorname{tr} P = p$。

証明のねらいは、$P$ が「射影らしい」性質を持つことの確認です。$P^2 = P$ は、一度射影した点は平面の上にあるので、もう一度射影しても動かない、という意味です。$P\boldsymbol{y}$ が射影であること自体は 定理 4.2 で示してあるので、ここでは性質を 1 つずつ計算で確かめます。

証明（命題 4.3）

1. $(X^{\mathsf{T}}X)^{-1}$ は対称行列の逆行列なので対称です（$(A^{-1})^{\mathsf{T}} = (A^{\mathsf{T}})^{-1}$ に $A^{\mathsf{T}}=A$ を代入）。よって $P^{\mathsf{T}} = X\bigl((X^{\mathsf{T}}X)^{-1}\bigr)^{\mathsf{T}}X^{\mathsf{T}} = P$ です。また

$$P^2 = X(X^{\mathsf{T}}X)^{-1}\underbrace{X^{\mathsf{T}}X(X^{\mathsf{T}}X)^{-1}}_{= I_p}X^{\mathsf{T}} = X(X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}} = P$$

です。

2. $\hat{\boldsymbol{w}} = (X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}\boldsymbol{y}$ は正規方程式の解ですから（両辺に $X^{\mathsf{T}}X$ を掛ければ確かめられます）、$P\boldsymbol{y} = X\hat{\boldsymbol{w}}$ であり、定理 4.2 の (2) よりこれは直交射影です。

3. $\hat{\boldsymbol{r}} = \boldsymbol{y}-P\boldsymbol{y} = (I_n-P)\boldsymbol{y}$ で、定理 4.2 の証明から $\hat{\boldsymbol{r}} \in (\operatorname{Im}X)^{\perp}$ です。$(I_n-P)^{\mathsf{T}} = I_n - P$、$(I_n-P)^2 = I_n - 2P + P^2 = I_n - P$ なのでこれも対称冪等であり、$\boldsymbol{u}\in(\operatorname{Im}X)^{\perp}$ に対しては 補題 4.1 より $X^{\mathsf{T}}\boldsymbol{u} = \boldsymbol{0}$、ゆえに $P\boldsymbol{u} = \boldsymbol{0}$、$(I_n-P)\boldsymbol{u} = \boldsymbol{u}$ となって、$(I_n-P)$ は $(\operatorname{Im}X)^{\perp}$ 上で恒等写像です。

4. トレースの巡回性 $\operatorname{tr}(AB) = \operatorname{tr}(BA)$ を $A = X$、$B = (X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}$ に使うと

$$\operatorname{tr}P = \operatorname{tr}\bigl((X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}X\bigr) = \operatorname{tr}I_p = p$$

です。

∎

$\hat{\boldsymbol{y}} = P\boldsymbol{y}$ が「$\boldsymbol{y}$ に帽子をかぶせる」ので、$P$ はハット行列と呼ばれます。統計学では $P$ の対角成分 $P_{ii}$ が「第 $i$ 観測が自分自身の予測値をどれだけ引っ張るか」を表す量（てこ比、leverage）として使われ、$\operatorname{tr}P = p$ は「$n$ 個の観測に $p$ 個分の自由度が費やされた」と読めます。

系 4.4（平方和の分解と決定係数）✓

$X$ の列に全成分 $1$ のベクトル $\boldsymbol{1}$ が含まれているとする。$\hat{\boldsymbol{y}} = X\hat{\boldsymbol{w}}$、$\hat{\boldsymbol{r}} = \boldsymbol{y}-\hat{\boldsymbol{y}}$、$\bar{y} = \frac{1}{n}\sum_i y_i$ とおくと、次が成り立つ。

1. $\sum_{i=1}^n \hat{r}_i = 0$、したがって $\frac{1}{n}\sum_i \hat{y}_i = \bar{y}$。
2. 平方和の分解

$$\underbrace{\sum_{i=1}^n (y_i-\bar{y})^2}_{\text{全変動 } S_{\mathrm{tot}}} = \underbrace{\sum_{i=1}^n (\hat{y}_i-\bar{y})^2}_{\text{回帰変動}} + \underbrace{\sum_{i=1}^n \hat{r}_i^{\,2}}_{\text{残差変動 } S_{\mathrm{res}}}$$

が成り立つ。とくに $S_{\mathrm{tot}} \ne 0$ のとき、決定係数 $R^2 = 1 - S_{\mathrm{res}}/S_{\mathrm{tot}}$ は $0 \le R^2 \le 1$ を満たす。

証明（系 4.4）

1. 定理 3.3 より $X^{\mathsf{T}}\hat{\boldsymbol{r}} = \boldsymbol{0}$ で、その $\boldsymbol{1}$ に対応する成分が $\langle \boldsymbol{1}, \hat{\boldsymbol{r}}\rangle = \sum_i \hat{r}_i = 0$ です。両辺を $n$ で割って $\bar{y} - \frac{1}{n}\sum_i\hat{y}_i = 0$ を得ます。

2. $\boldsymbol{y}-\bar{y}\boldsymbol{1} = (\hat{\boldsymbol{y}}-\bar{y}\boldsymbol{1}) + \hat{\boldsymbol{r}}$ と分解します。$\boldsymbol{1} \in \operatorname{Im}X$ かつ $\hat{\boldsymbol{y}}\in\operatorname{Im}X$ なので $\hat{\boldsymbol{y}}-\bar{y}\boldsymbol{1}\in\operatorname{Im}X$ であり、定理 4.2 より $\hat{\boldsymbol{r}} \in (\operatorname{Im}X)^{\perp}$ ですから、この 2 つは直交します。ピタゴラスの定理より主張の等式が出ます。$R^2 = (\text{回帰変動})/S_{\mathrm{tot}}$ であり、両変動とも非負で和が $S_{\mathrm{tot}}$ なので $0\le R^2\le 1$ です。

∎

例 4.5（3 点への直線当てはめを射影として見る）

$n = 3$、$d = 1$ で、データを $(x_i, y_i) = (0,1), (1,1), (2,4)$ とします。計画行列と観測ベクトルは

$$X = \begin{pmatrix} 1 & 0 \\ 1 & 1 \\ 1 & 2\end{pmatrix},\qquad \boldsymbol{y} = \begin{pmatrix} 1 \\ 1 \\ 4\end{pmatrix}$$

です。$\operatorname{Im}X$ は $\mathbb{R}^3$ の中の 2 次元平面で、$\boldsymbol{y}$ はその上にありません（3 点は一直線上にないため）。正規方程式の材料を計算すると

$$X^{\mathsf{T}}X = \begin{pmatrix} 3 & 3 \\ 3 & 5\end{pmatrix},\qquad X^{\mathsf{T}}\boldsymbol{y} = \begin{pmatrix} 1+1+4 \\ 0\cdot 1 + 1\cdot 1 + 2\cdot 4\end{pmatrix} = \begin{pmatrix} 6 \\ 9\end{pmatrix}$$

です。$\det(X^{\mathsf{T}}X) = 15-9 = 6 \ne 0$ なので逆行列が存在し、

$$\hat{\boldsymbol{w}} = \frac{1}{6}\begin{pmatrix} 5 & -3 \\ -3 & 3\end{pmatrix}\begin{pmatrix} 6 \\ 9\end{pmatrix} = \frac{1}{6}\begin{pmatrix} 30-27 \\ -18+27\end{pmatrix} = \begin{pmatrix} 0.5 \\ 1.5\end{pmatrix}$$

すなわち $\hat{y} = 0.5 + 1.5x$ を得ます。予測と残差は

$$\hat{\boldsymbol{y}} = X\hat{\boldsymbol{w}} = (0.5,\ 2.0,\ 3.5)^{\mathsf{T}},\qquad \hat{\boldsymbol{r}} = \boldsymbol{y}-\hat{\boldsymbol{y}} = (0.5,\ -1.0,\ 0.5)^{\mathsf{T}}$$

です。直交性を確かめます。$\langle \boldsymbol{1}, \hat{\boldsymbol{r}}\rangle = 0.5-1.0+0.5 = 0$、$\langle (0,1,2)^{\mathsf{T}}, \hat{\boldsymbol{r}}\rangle = 0\cdot 0.5 + 1\cdot(-1.0) + 2\cdot 0.5 = 0$ となり、確かに残差は $X$ の両方の列と直交しています。ピタゴラスの定理も

$$\lVert\boldsymbol{y}\rVert^2 = 1+1+16 = 18,\quad \lVert\hat{\boldsymbol{y}}\rVert^2 = 0.25+4+12.25 = 16.5,\quad \lVert\hat{\boldsymbol{r}}\rVert^2 = 0.25+1+0.25 = 1.5$$

より $18 = 16.5 + 1.5$ と成立しています。ハット行列も具体的に書けて

$$P = X(X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}} = \frac{1}{6}\begin{pmatrix} 5 & 2 & -1 \\ 2 & 2 & 2 \\ -1 & 2 & 5\end{pmatrix}$$

となります。実際 $P\boldsymbol{y} = \frac{1}{6}(5+2-4,\ 2+2+8,\ -1+2+20)^{\mathsf{T}} = (0.5, 2.0, 3.5)^{\mathsf{T}} = \hat{\boldsymbol{y}}$ です。対称性は見てのとおりで、$\operatorname{tr}P = (5+2+5)/6 = 2 = p$ も 命題 4.3 の 4 と合致します。

5. 解はいつ一意になるか：ランク条件

定理 4.2 により、予測 $\hat{\boldsymbol{y}} = X\hat{\boldsymbol{w}}$ はつねに一意でした。しかし係数 $\hat{\boldsymbol{w}}$ そのものが一意かどうかは別問題です。定理 3.3 の最後の主張から、解集合は $\hat{\boldsymbol{w}} + \operatorname{Ker}X$ ですから、一意性は $\operatorname{Ker}X = \{\boldsymbol{0}\}$ と同値です。これを $X^{\mathsf{T}}X$ の言葉に翻訳します。

補題 5.1（グラム行列の核）✓

任意の $X\in\mathbb{R}^{n\times p}$ に対し

$$\operatorname{Ker}(X^{\mathsf{T}}X) = \operatorname{Ker}X, \qquad \operatorname{rank}(X^{\mathsf{T}}X) = \operatorname{rank}X$$

が成り立つ。

証明（補題 5.1）

$X\boldsymbol{v} = \boldsymbol{0}$ ならば両辺に左から $X^{\mathsf{T}}$ を掛けて $X^{\mathsf{T}}X\boldsymbol{v} = \boldsymbol{0}$ ですから $\operatorname{Ker}X \subseteq \operatorname{Ker}(X^{\mathsf{T}}X)$ です。逆に $X^{\mathsf{T}}X\boldsymbol{v} = \boldsymbol{0}$ とすると、左から $\boldsymbol{v}^{\mathsf{T}}$ を掛けて

$$0 = \boldsymbol{v}^{\mathsf{T}}X^{\mathsf{T}}X\boldsymbol{v} = (X\boldsymbol{v})^{\mathsf{T}}(X\boldsymbol{v}) = \lVert X\boldsymbol{v}\rVert^2$$

となります。ノルムが $0$ になるのはゼロベクトルだけなので $X\boldsymbol{v} = \boldsymbol{0}$、すなわち $\operatorname{Ker}(X^{\mathsf{T}}X)\subseteq\operatorname{Ker}X$ です。よって両者は一致します。階数については、$X$ も $X^{\mathsf{T}}X$ も列数が $p$ なので、次元定理

$$\operatorname{rank}A = p - \dim\operatorname{Ker}A$$

を両者に適用すれば、核の次元が等しいことから階数も等しくなります。

∎

系 5.2（最小二乗解の一意性）✓

$X\in\mathbb{R}^{n\times p}$、$\boldsymbol{y}\in\mathbb{R}^n$ とする。次の 3 条件は同値である。

1. $X$ の列 $p$ 本は一次独立、すなわち $\operatorname{rank}X = p$。
2. $X^{\mathsf{T}}X$ は正則。
3. 最小二乗解がただ 1 つ存在する。

このとき最小二乗解は

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

で与えられる。さらに、これらが成り立つためには $n \ge p$ が必要である。

証明（系 5.2）

1 $\Leftrightarrow$ 2。 $X^{\mathsf{T}}X$ は $p$ 次正方行列なので、正則であることは $\operatorname{rank}(X^{\mathsf{T}}X) = p$ と同値です。補題 5.1 よりこれは $\operatorname{rank}X = p$ と同値で、これは列が一次独立であることにほかなりません。

1 $\Leftrightarrow$ 3。 定理 4.2 より最小二乗解は必ず存在し、定理 3.3 より解集合は $\hat{\boldsymbol{w}} + \operatorname{Ker}X$ です。したがって解が 1 つであることは $\operatorname{Ker}X = \{\boldsymbol{0}\}$ と同値で、次元定理よりこれは $\operatorname{rank}X = p$ と同値です。

解の公式。 2 のもとで $\boldsymbol{w} = (X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}\boldsymbol{y}$ とおくと $X^{\mathsf{T}}X\boldsymbol{w} = X^{\mathsf{T}}\boldsymbol{y}$ が成り立つので、これは正規方程式の解であり、定理 3.3 より最小二乗解です。一意性はいま示したとおりです。

$n\ge p$ の必要性。 $\operatorname{rank}X \le \min(n,p)$ ですから、$\operatorname{rank}X = p$ なら $p \le n$ です。

∎

データ数が特徴量数より少ない（$n < p$）場合には、この意味で解は決して一意になりません。深層学習で日常的に起きる「パラメータのほうが多い」状況が、すでに線形回帰の段階で顔を出しているわけです。

例 5.3（特徴量が重複しているとき（多重共線性））

例 4.5 と同じデータ $(x_i, y_i) = (0,1), (1,1), (2,4)$ に対し、「$x$ の 2 倍」という無意味な特徴量 $x' = 2x$ を追加してみます。計画行列は

$$X = \begin{pmatrix} 1 & 0 & 0 \\ 1 & 1 & 2 \\ 1 & 2 & 4\end{pmatrix}$$

となり、第 3 列は第 2 列のちょうど 2 倍なので $\operatorname{rank}X = 2 < 3 = p$ です。核は

$$\operatorname{Ker}X = \{\, t(0, 2, -1)^{\mathsf{T}} : t\in\mathbb{R}\,\}$$

です（実際 $0\cdot\boldsymbol{1} + 2\boldsymbol{x} - 1\cdot(2\boldsymbol{x}) = \boldsymbol{0}$）。例 4.5 で求めた $(0.5, 1.5)$ に第 3 座標 $0$ を足した $\hat{\boldsymbol{w}}_0 = (0.5, 1.5, 0)^{\mathsf{T}}$ は最小二乗解ですから、定理 3.3 より解集合は

$$\{\,(0.5,\ 1.5 - 2t,\ t)^{\mathsf{T}} : t\in\mathbb{R}\,\}$$

という直線全体になります。たとえば $t = 0.75$ は $(0.5, 0, 0.75)$、$t = -1$ は $(0.5, 3.5, -1)$ で、係数の見た目はまったく違います。「$x$ の効果は $1.5$」なのか「$x$ の効果は $0$ で $x'$ の効果が $0.75$」なのかは、データからは決められません。それでも予測は

$$X(0.5,\ 1.5-2t,\ t)^{\mathsf{T}} = (0.5,\ 0.5 + (1.5-2t) + 2t,\ 0.5 + 2(1.5-2t)+4t)^{\mathsf{T}} = (0.5,\ 2,\ 3.5)^{\mathsf{T}}$$

と $t$ によらず一定で、定理 4.2 の (2) のとおりです。

注意

実務で厄介なのは、厳密に一次従属な場合よりほぼ一次従属な場合です。$\operatorname{rank}X = p$ ではあるが $X^{\mathsf{T}}X$ の最小固有値がほぼ $0$ という状況では、$\hat{\boldsymbol{w}}$ は一意に決まるものの、$\boldsymbol{y}$ のわずかな変化で激しく振れます。身長を cm と inch の両方で特徴量に入れる、カテゴリ変数をダミー化するときに全水準の列と切片列を同時に入れる（いわゆるダミー変数の罠）といった操作は、この状態か、その極限であるランク落ちを作ります（ダミー変数の罠は後者で、全水準のダミー列の和がちょうど $\boldsymbol{1}$ になるため厳密に一次従属です）。対処は、冗長な列を落とすか、演習 8.3 のリッジ回帰のように正則化することです。

6. 手を動かして解く

6.1. 単回帰の閉じた公式

$d = 1$（説明変数 1 個）の場合を、一般式まで書き下します。$p = 2$ で

$$X^{\mathsf{T}}X = \begin{pmatrix} n & \sum_i x_i \\ \sum_i x_i & \sum_i x_i^2\end{pmatrix}, \qquad X^{\mathsf{T}}\boldsymbol{y} = \begin{pmatrix} \sum_i y_i \\ \sum_i x_i y_i\end{pmatrix}$$

です。$\bar{x} = \frac1n\sum_i x_i$、$\bar{y} = \frac1n\sum_i y_i$、$S_{xx} = \sum_i (x_i-\bar{x})^2$、$S_{xy} = \sum_i (x_i-\bar{x})(y_i-\bar{y})$ とおくと、展開して

$$\det(X^{\mathsf{T}}X) = n\sum_i x_i^2 - \Bigl(\sum_i x_i\Bigr)^2 = n\Bigl(\sum_i x_i^2 - n\bar{x}^2\Bigr) = nS_{xx}$$

なので、$S_{xx}\ne 0$（$x_i$ が全部同じ値ではない）なら 系 5.2 が使えて

$$\hat{\boldsymbol{w}} = \frac{1}{nS_{xx}}\begin{pmatrix} \sum_i x_i^2 & -\sum_i x_i \\ -\sum_i x_i & n\end{pmatrix}\begin{pmatrix} n\bar{y} \\ \sum_i x_i y_i\end{pmatrix}$$

となります。第 2 成分は $\bigl(n\sum_i x_iy_i - n\bar{y}\sum_i x_i\bigr)/(nS_{xx}) = (\sum_i x_iy_i - n\bar{x}\bar{y})/S_{xx} = S_{xy}/S_{xx}$ です。第 1 成分も同様に計算すると $(\bar{y}\sum_i x_i^2 - \bar{x}\sum_i x_iy_i)/S_{xx}$ となり、これは $\bar{y}S_{xx} - \bar{x}S_{xy} = \bar{y}\sum_i x_i^2 - n\bar{x}^2\bar{y} - \bar{x}\sum_i x_iy_i + n\bar{x}^2\bar{y} = \bar{y}\sum_i x_i^2 - \bar{x}\sum_i x_iy_i$ に一致します。まとめると、高校でも見る形

$$\hat{w}_1 = \frac{S_{xy}}{S_{xx}}, \qquad \hat{b} = \bar{y} - \hat{w}_1\bar{x}$$

が得られます。第 2 式は「回帰直線は必ず重心 $(\bar{x},\bar{y})$ を通る」と読めます。これは 系 4.4 の 1（残差の和が $0$）の言い換えでもあります。

例 6.1（5 点の単回帰を最後まで計算する）

データを $(x_i,y_i) = (1,2), (2,3), (3,5), (4,4), (5,6)$ とします。$n = 5$ で

$$\sum_i x_i = 15,\quad \sum_i y_i = 20,\quad \sum_i x_i^2 = 55,\quad \sum_i x_iy_i = 2+6+15+16+30 = 69$$

ですから $\bar{x} = 3$、$\bar{y} = 4$ です。正規方程式は

$$\begin{pmatrix} 5 & 15 \\ 15 & 55\end{pmatrix}\begin{pmatrix} b \\ w_1\end{pmatrix} = \begin{pmatrix} 20 \\ 69\end{pmatrix}$$

です。$\det = 5\cdot 55 - 15^2 = 275-225 = 50 \ne 0$ なので

$$\begin{pmatrix} \hat{b} \\ \hat{w}_1\end{pmatrix} = \frac{1}{50}\begin{pmatrix} 55 & -15 \\ -15 & 5\end{pmatrix}\begin{pmatrix} 20 \\ 69\end{pmatrix} = \frac{1}{50}\begin{pmatrix} 1100 - 1035 \\ -300 + 345\end{pmatrix} = \frac{1}{50}\begin{pmatrix} 65 \\ 45\end{pmatrix} = \begin{pmatrix} 1.3 \\ 0.9\end{pmatrix}$$

すなわち $\hat{y} = 1.3 + 0.9x$ です。公式でも確かめます。$S_{xx} = 4+1+0+1+4 = 10$、$S_{xy} = (-2)(-2)+(-1)(-1)+0\cdot 1+1\cdot 0+2\cdot 2 = 9$ なので $\hat{w}_1 = 9/10 = 0.9$、$\hat{b} = 4 - 0.9\cdot 3 = 1.3$ と一致しました。

予測値は $\hat{\boldsymbol{y}} = (2.2,\ 3.1,\ 4.0,\ 4.9,\ 5.8)^{\mathsf{T}}$、残差は

$$\hat{\boldsymbol{r}} = (-0.2,\ -0.1,\ 1.0,\ -0.9,\ 0.2)^{\mathsf{T}}$$

です。直交条件を検算すると $\sum_i \hat{r}_i = -0.2-0.1+1.0-0.9+0.2 = 0$、$\sum_i x_i\hat{r}_i = -0.2-0.2+3.0-3.6+1.0 = 0$ でどちらも $0$ になります。平方和は $S_{\mathrm{res}} = 0.04+0.01+1.00+0.81+0.04 = 1.90$、$S_{\mathrm{tot}} = 4+1+1+0+4 = 10$ なので、系 4.4 より回帰変動は $10-1.9 = 8.1$ で、これは $S_{xy}^2/S_{xx} = 81/10 = 8.1$ と一致します。決定係数は $R^2 = 1 - 1.9/10 = 0.81$ です。

5 点の単回帰。直線 y = 1.3 + 0.9x は重心 (3, 4) を通り、点線の縦線が残差である。上に出る残差と下に出る残差が打ち消し合って、残差の和は 0 になる。

同じ計算を NumPy で書くと次のようになります。逆行列を作らずに連立一次方程式として解いている点に注目してください。

import numpy as np

x = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
y = np.array([2.0, 3.0, 5.0, 4.0, 6.0])
X = np.column_stack([np.ones_like(x), x])        # 計画行列 (5, 2)

# 正規方程式 X^T X w = X^T y を直接解く
w_ne = np.linalg.solve(X.T @ X, X.T @ y)         # [1.3, 0.9]

# 数値的にはこちらが推奨（内部で特異値分解を使う）
w_ls, *_ = np.linalg.lstsq(X, y, rcond=None)     # [1.3, 0.9]

r = y - X @ w_ls
print(X.T @ r)          # [0. 0.] にごく近い値（丸め誤差の範囲）
print(1 - r @ r / ((y - y.mean()) ** 2).sum())   # 0.81

6.2. 多項式当てはめ

定義 2.2 のとおり、特徴写像を取り替えるだけで曲線の当てはめになります。道具立ては何も変わりません。

例 6.2（2 次多項式を当てはめる）

データを $(x_i, y_i) = (-2,6), (-1,2), (0,2), (1,2), (2,5)$ とし、$\varphi(x) = (1, x, x^2)^{\mathsf{T}}$、すなわち $\hat{y} = w_0 + w_1x + w_2x^2$ を当てはめます。$p = 3$ で、必要な和は

$$\sum_i x_i = 0,\quad \sum_i x_i^2 = 10,\quad \sum_i x_i^3 = 0,\quad \sum_i x_i^4 = 34$$

（$x_i$ が $0$ について対称なので奇数乗の和が消えます）、および

$$\sum_i y_i = 17,\quad \sum_i x_iy_i = -12-2+0+2+10 = -2,\quad \sum_i x_i^2y_i = 24+2+0+2+20 = 48$$

です。したがって正規方程式は

$$\begin{pmatrix} 5 & 0 & 10 \\ 0 & 10 & 0 \\ 10 & 0 & 34 \end{pmatrix} \begin{pmatrix} w_0 \\ w_1 \\ w_2\end{pmatrix} = \begin{pmatrix} 17 \\ -2 \\ 48\end{pmatrix}$$

となります。第 2 行はただちに $\hat{w}_1 = -0.2$ を与えます。第 1 行と第 3 行は $w_0, w_2$ だけの連立方程式で、第 1 行を 2 倍して第 3 行から引くと

$$(34 - 20)w_2 = 48 - 34, \qquad \text{すなわち}\quad 14w_2 = 14$$

なので $\hat{w}_2 = 1$、そして第 1 行から $\hat{w}_0 = (17-10)/5 = 1.4$ です。当てはめた曲線は

$$\hat{y} = 1.4 - 0.2x + x^2$$

です。予測値は $\hat{\boldsymbol{y}} = (5.8,\ 2.6,\ 1.4,\ 2.2,\ 5.0)^{\mathsf{T}}$、残差は $\hat{\boldsymbol{r}} = (0.2,\ -0.6,\ 0.6,\ -0.2,\ 0)^{\mathsf{T}}$ です。ここでも直交条件を検算すると、3 本すべてが成り立ちます。

$$\begin{aligned} \sum_i \hat{r}_i &= 0.2-0.6+0.6-0.2+0 = 0,\\ \sum_i x_i\hat{r}_i &= -0.4+0.6+0-0.2+0 = 0,\\ \sum_i x_i^2\hat{r}_i &= 0.8-0.6+0-0.2+0 = 0. \end{aligned}$$

残差平方和は $S_{\mathrm{res}} = 0.04+0.36+0.36+0.04+0 = 0.8$、全変動は $\bar{y} = 3.4$ より $S_{\mathrm{tot}} = 6.76+1.96+1.96+1.96+2.56 = 15.2$ なので $R^2 = 1-0.8/15.2 \approx 0.947$ です。曲線を当てはめたのに、解いたのは 3 元連立一次方程式にすぎません。これが「パラメータについて線形」であることの威力です。

7. 数値計算上の注意

7.1. 解法の全体像

ここまでの結果を 1 枚にまとめておきます。

最小二乗問題から解までの流れ。分岐はランク条件で、予測ベクトルはどちらの枝でも一意に定まる。

$\operatorname{rank}X = p$ のとき $X^{\mathsf{T}}X$ は対称かつ正定値です。実際 補題 5.1 より $\boldsymbol{v}\ne\boldsymbol{0}$ なら $X\boldsymbol{v}\ne\boldsymbol{0}$ なので $\boldsymbol{v}^{\mathsf{T}}X^{\mathsf{T}}X\boldsymbol{v} = \lVert X\boldsymbol{v}\rVert^2 > 0$ です。正定値対称行列にはコレスキー分解（行列を下三角行列とその転置の積に分ける方法）$X^{\mathsf{T}}X = LL^{\mathsf{T}}$（$L$ は下三角）が存在するので、実装は「$X^{\mathsf{T}}X$ と $X^{\mathsf{T}}\boldsymbol{y}$ を作る、コレスキー分解する、三角方程式を 2 回解く」で終わります。

ヒント

np.linalg.inv(X.T @ X) @ X.T @ y と書きたくなりますが、やめてください。逆行列を陽に作るのは計算量の面でも精度の面でも損です。連立一次方程式として解く（np.linalg.solve）か、次に述べる QR・特異値分解を使う（np.linalg.lstsq）のが定石です。

7.2. 条件数は二乗される

この節で言いたいのは、$X$ がもつ「ゆがみ具合」（最も大きく伸ばす向きと最も小さく伸ばす向きの比）が、$X^{\mathsf{T}}X$ をつくると二乗になってしまう、ということです。特異値は行列が各方向をどれだけ伸ばすかを表す数、特異値分解はそれを取り出す分解だと思ってください。

正規方程式には理論上きれいでも数値計算上は不利な点があります。行列 $A$ の条件数を $\kappa_2(A) = \sigma_{\max}(A)/\sigma_{\min}(A)$（特異値の比）で測ると、フルランクの $X$ に対して

$$\kappa_2(X^{\mathsf{T}}X) = \kappa_2(X)^2$$

が成り立ちます。理由は、$X$ の特異値分解 $X = U\Sigma V^{\mathsf{T}}$ を使うと $X^{\mathsf{T}}X = V\Sigma^{\mathsf{T}}\Sigma V^{\mathsf{T}}$ となり、$X^{\mathsf{T}}X$ の固有値が $X$ の特異値の二乗になるからです（対称行列の対角化については スペクトル定理 の 実対称行列の直交対角化（系 4.3）［スペクトル定理］ を参照してください）。条件数は「入力の相対誤差が解の相対誤差に何倍に増幅されるか」の目安ですから、$X^{\mathsf{T}}X$ を作った瞬間に増幅率が二乗されることになります。

例 7.1（中心化するだけで条件数が改善する）

例 6.1 のデータでは $X^{\mathsf{T}}X = \begin{pmatrix} 5 & 15 \\ 15 & 55\end{pmatrix}$ でした。対称 $2\times 2$ 行列なので固有値は特性方程式 $\lambda^2 - 60\lambda + 50 = 0$（トレース $60$、行列式 $50$）から

$$\lambda = 30 \pm \sqrt{900-50} = 30 \pm \sqrt{850}$$

すなわち $\lambda_{\max} \approx 59.155$、$\lambda_{\min} \approx 0.845$ です。よって $\kappa_2(X^{\mathsf{T}}X) \approx 70.0$、$\kappa_2(X) \approx \sqrt{70} \approx 8.4$ です。

ここで $x$ を中心化して第 2 列を $x_i - \bar{x} = -2,-1,0,1,2$ に置き換えると、第 1 列との内積が $\sum_i (x_i-\bar{x}) = 0$ になるので

$$X_{\mathrm{c}}^{\mathsf{T}}X_{\mathrm{c}} = \begin{pmatrix} 5 & 0 \\ 0 & 10\end{pmatrix}$$

と対角になり、$\kappa_2(X_{\mathrm{c}}^{\mathsf{T}}X_{\mathrm{c}}) = 10/5 = 2$ まで下がります。データの中身は同じで、座標の取り方を変えただけです。特徴量の中心化とスケーリングが前処理として推奨される理由の一つがこれです（中心化しても当てはめた直線そのものは変わりません。演習 8.2 を見てください）。

注意 7.2（QR 分解による解法）

$\operatorname{rank}X = p$ のとき、$X$ は $X = QR$（$Q\in\mathbb{R}^{n\times p}$ は $Q^{\mathsf{T}}Q = I_p$、$R\in\mathbb{R}^{p\times p}$ は対角成分が正の上三角行列）と分解できます。これは グラム・シュミットの直交化（定理 6.1）［内積空間とグラム・シュミット直交化］ そのものです。これを正規方程式に代入すると

$$X^{\mathsf{T}}X = R^{\mathsf{T}}Q^{\mathsf{T}}QR = R^{\mathsf{T}}R, \qquad X^{\mathsf{T}}\boldsymbol{y} = R^{\mathsf{T}}Q^{\mathsf{T}}\boldsymbol{y}$$

なので、正規方程式は $R^{\mathsf{T}}R\boldsymbol{w} = R^{\mathsf{T}}Q^{\mathsf{T}}\boldsymbol{y}$ となります。$R$ は正則（対角成分が正）なので $R^{\mathsf{T}}$ も正則で、両辺に $(R^{\mathsf{T}})^{-1}$ を掛ければ

$$R\hat{\boldsymbol{w}} = Q^{\mathsf{T}}\boldsymbol{y}$$

という上三角の連立方程式に化けます。後退代入で解けて、しかも $X^{\mathsf{T}}X$ を作らないので条件数は $\kappa_2(X)$ のままです。ハット行列も $P = QQ^{\mathsf{T}}$ と簡単になります。ランク落ちの可能性がある場合には、さらに頑健な特異値分解による解法（本記事末尾の Appendix）が使われます。

7.3. この先へ

閉じた解が使えなくなる典型的な場面と、その先の道具を挙げておきます。

状況	起きること	対処
$p$ が数万以上	$p\times p$ 行列の分解が $O(p^3)$ で重い	勾配降下法・共役勾配法
$\operatorname{rank}X < p$、または列がほぼ従属	解が不定、係数が暴れる	リッジ回帰（演習 8.3）、擬似逆行列（Appendix）
目的変数が $0$ か $1$	二乗誤差が不自然、確率にならない	ロジスティック回帰
特徴量そのものを作りたい	何を $\varphi$ に選ぶかが問題になる	主成分分析、ニューラルネットワークと逆伝播
係数の不確かさを知りたい	点推定だけでは足りない	確率論とベイズ統計の役割

なお統計学の側からは、誤差 $\varepsilon_i$ が平均 $0$・分散 $\sigma^2$ で無相関という仮定のもとで、最小二乗推定量が「不偏な線形推定量のなかで分散最小」であること（ガウス・マルコフの定理）が示されます。本記事の枠組みでは確率を仮定していないので触れませんが、証明はここで作ったハット行列と直交性だけで書けます。

8. 演習

演習 8.1易

計画行列 $X$ の列に $\boldsymbol{1} = (1,\dots,1)^{\mathsf{T}}$ が含まれるならば、最小二乗解の残差は $\sum_i \hat{r}_i = 0$ を満たすことを示してください。また、切片を持たないモデル $\hat{y} = w x$（計画行列は $\boldsymbol{x}$ の 1 列のみ）ではこれが成り立たないことを、具体的なデータで確かめてください。

解答

前半。定理 3.3 より最小二乗解は $X^{\mathsf{T}}(\boldsymbol{y}-X\hat{\boldsymbol{w}}) = X^{\mathsf{T}}\hat{\boldsymbol{r}} = \boldsymbol{0}$ を満たします。$X^{\mathsf{T}}\hat{\boldsymbol{r}}$ の第 $j$ 成分は $X$ の第 $j$ 列と $\hat{\boldsymbol{r}}$ の内積です。第 $j$ 列が $\boldsymbol{1}$ ならその成分は $\langle\boldsymbol{1},\hat{\boldsymbol{r}}\rangle = \sum_i \hat{r}_i$ であり、これが $0$ になります。

後半。データを $(x_1,y_1) = (1,1)$、$(x_2,y_2) = (2,0)$ とします。計画行列は $X = (1, 2)^{\mathsf{T}}$（$n=2$、$p=1$）で、$X^{\mathsf{T}}X = 1^2+2^2 = 5$、$X^{\mathsf{T}}\boldsymbol{y} = 1\cdot 1 + 2\cdot 0 = 1$ なので $\hat{w} = 1/5 = 0.2$ です。予測は $\hat{\boldsymbol{y}} = (0.2, 0.4)^{\mathsf{T}}$、残差は $\hat{\boldsymbol{r}} = (0.8, -0.4)^{\mathsf{T}}$ で、$\sum_i \hat{r}_i = 0.4 \ne 0$ です。もちろん $\langle \boldsymbol{x}, \hat{\boldsymbol{r}}\rangle = 1\cdot 0.8 + 2\cdot(-0.4) = 0$ は成り立っています。直交するのは「$X$ の列」であって、モデルに含めていない $\boldsymbol{1}$ とは直交しません。

演習 8.2標準

$X = (\boldsymbol{1} \ \ \tilde{X})$、$\tilde{X}\in\mathbb{R}^{n\times d}$ と分割し、パラメータも $\boldsymbol{w} = (b, \boldsymbol{u})$（$\boldsymbol{u}\in\mathbb{R}^d$）と分ける。列平均を $\bar{\boldsymbol{x}} = \frac1n\tilde{X}^{\mathsf{T}}\boldsymbol{1}\in\mathbb{R}^d$、中心化した行列とベクトルを $X_{\mathrm{c}} = \tilde{X}-\boldsymbol{1}\bar{\boldsymbol{x}}^{\mathsf{T}}$、$\boldsymbol{y}_{\mathrm{c}} = \boldsymbol{y}-\bar{y}\boldsymbol{1}$ とおく。このとき正規方程式が

$$\hat{b} = \bar{y}-\bar{\boldsymbol{x}}^{\mathsf{T}}\hat{\boldsymbol{u}}, \qquad X_{\mathrm{c}}^{\mathsf{T}}X_{\mathrm{c}}\,\hat{\boldsymbol{u}} = X_{\mathrm{c}}^{\mathsf{T}}\boldsymbol{y}_{\mathrm{c}}$$

と同値であることを示してください。

解答

正規方程式 $X^{\mathsf{T}}X\boldsymbol{w} = X^{\mathsf{T}}\boldsymbol{y}$ をブロックで書きます。$\boldsymbol{1}^{\mathsf{T}}\boldsymbol{1} = n$、$\boldsymbol{1}^{\mathsf{T}}\tilde{X} = n\bar{\boldsymbol{x}}^{\mathsf{T}}$、$\boldsymbol{1}^{\mathsf{T}}\boldsymbol{y} = n\bar{y}$ なので

$$\begin{pmatrix} n & n\bar{\boldsymbol{x}}^{\mathsf{T}} \\ n\bar{\boldsymbol{x}} & \tilde{X}^{\mathsf{T}}\tilde{X}\end{pmatrix} \begin{pmatrix} b \\ \boldsymbol{u}\end{pmatrix} = \begin{pmatrix} n\bar{y} \\ \tilde{X}^{\mathsf{T}}\boldsymbol{y}\end{pmatrix}$$

です。第 1 ブロック行は $nb + n\bar{\boldsymbol{x}}^{\mathsf{T}}\boldsymbol{u} = n\bar{y}$、すなわち $b = \bar{y}-\bar{\boldsymbol{x}}^{\mathsf{T}}\boldsymbol{u}$ です。これを第 2 ブロック行 $n\bar{\boldsymbol{x}}b + \tilde{X}^{\mathsf{T}}\tilde{X}\boldsymbol{u} = \tilde{X}^{\mathsf{T}}\boldsymbol{y}$ に代入すると

$$\bigl(\tilde{X}^{\mathsf{T}}\tilde{X} - n\bar{\boldsymbol{x}}\bar{\boldsymbol{x}}^{\mathsf{T}}\bigr)\boldsymbol{u} = \tilde{X}^{\mathsf{T}}\boldsymbol{y} - n\bar{\boldsymbol{x}}\bar{y}$$

を得ます。あとは左辺と右辺の係数が中心化された量であることを確かめます。$\boldsymbol{1}^{\mathsf{T}}\boldsymbol{1} = n$ と $\boldsymbol{1}^{\mathsf{T}}\tilde{X} = n\bar{\boldsymbol{x}}^{\mathsf{T}}$ より

$$\begin{aligned} X_{\mathrm{c}}^{\mathsf{T}}X_{\mathrm{c}} &= (\tilde{X}-\boldsymbol{1}\bar{\boldsymbol{x}}^{\mathsf{T}})^{\mathsf{T}}(\tilde{X}-\boldsymbol{1}\bar{\boldsymbol{x}}^{\mathsf{T}}) \\ &= \tilde{X}^{\mathsf{T}}\tilde{X} - \tilde{X}^{\mathsf{T}}\boldsymbol{1}\bar{\boldsymbol{x}}^{\mathsf{T}} - \bar{\boldsymbol{x}}\boldsymbol{1}^{\mathsf{T}}\tilde{X} + \bar{\boldsymbol{x}}(\boldsymbol{1}^{\mathsf{T}}\boldsymbol{1})\bar{\boldsymbol{x}}^{\mathsf{T}} \\ &= \tilde{X}^{\mathsf{T}}\tilde{X} - n\bar{\boldsymbol{x}}\bar{\boldsymbol{x}}^{\mathsf{T}} - n\bar{\boldsymbol{x}}\bar{\boldsymbol{x}}^{\mathsf{T}} + n\bar{\boldsymbol{x}}\bar{\boldsymbol{x}}^{\mathsf{T}} = \tilde{X}^{\mathsf{T}}\tilde{X} - n\bar{\boldsymbol{x}}\bar{\boldsymbol{x}}^{\mathsf{T}} \end{aligned}$$

であり、同様に

$$X_{\mathrm{c}}^{\mathsf{T}}\boldsymbol{y}_{\mathrm{c}} = \tilde{X}^{\mathsf{T}}\boldsymbol{y} - \bar{y}\tilde{X}^{\mathsf{T}}\boldsymbol{1} - \bar{\boldsymbol{x}}\boldsymbol{1}^{\mathsf{T}}\boldsymbol{y} + \bar{y}\bar{\boldsymbol{x}}\boldsymbol{1}^{\mathsf{T}}\boldsymbol{1} = \tilde{X}^{\mathsf{T}}\boldsymbol{y} - n\bar{y}\bar{\boldsymbol{x}}$$

です。よって主張の 2 式が得られ、逆にこの 2 式から元のブロック方程式が復元できるので同値です。つまり、まずデータを中心化して $d$ 元の方程式を解き、最後に切片を $\hat{b} = \bar{y}-\bar{\boldsymbol{x}}^{\mathsf{T}}\hat{\boldsymbol{u}}$ で求めればよいことになります。$d = 1$ のときこれは 例 6.1 の公式にほかなりません。

演習 8.3標準

$\lambda > 0$ とし、リッジ回帰の目的関数

$$L_{\lambda}(\boldsymbol{w}) = \lVert \boldsymbol{y}-X\boldsymbol{w}\rVert^2 + \lambda\lVert\boldsymbol{w}\rVert^2$$

を考える。(1) $X$ のランクによらず $X^{\mathsf{T}}X+\lambda I_p$ が正則であることを示してください。(2) $L_{\lambda}$ の最小点がただ 1 つ存在し、$\hat{\boldsymbol{w}}_{\lambda} = (X^{\mathsf{T}}X+\lambda I_p)^{-1}X^{\mathsf{T}}\boldsymbol{y}$ で与えられることを示してください。

解答

(1)。 $\boldsymbol{v}\ne\boldsymbol{0}$ に対し

$$\boldsymbol{v}^{\mathsf{T}}(X^{\mathsf{T}}X+\lambda I_p)\boldsymbol{v} = \lVert X\boldsymbol{v}\rVert^2 + \lambda\lVert\boldsymbol{v}\rVert^2 \ge \lambda\lVert\boldsymbol{v}\rVert^2 > 0$$

です。もし $(X^{\mathsf{T}}X+\lambda I_p)\boldsymbol{v} = \boldsymbol{0}$ となる $\boldsymbol{v}\ne\boldsymbol{0}$ があれば左辺が $0$ になって矛盾するので、核は $\{\boldsymbol{0}\}$、すなわち正則です。

(2)。 拡大した計画行列と観測ベクトル

$$X_{\lambda} = \begin{pmatrix} X \\ \sqrt{\lambda}\,I_p\end{pmatrix}\in\mathbb{R}^{(n+p)\times p}, \qquad \boldsymbol{y}_{\lambda} = \begin{pmatrix} \boldsymbol{y} \\ \boldsymbol{0}\end{pmatrix}\in\mathbb{R}^{n+p}$$

を考えます。ブロックごとにノルムを計算すると

$$\lVert \boldsymbol{y}_{\lambda}-X_{\lambda}\boldsymbol{w}\rVert^2 = \lVert\boldsymbol{y}-X\boldsymbol{w}\rVert^2 + \lVert\sqrt{\lambda}\boldsymbol{w}\rVert^2 = L_{\lambda}(\boldsymbol{w})$$

なので、リッジ回帰は拡大データに対する通常の最小二乗問題です。$X_{\lambda}^{\mathsf{T}}X_{\lambda} = X^{\mathsf{T}}X+\lambda I_p$、$X_{\lambda}^{\mathsf{T}}\boldsymbol{y}_{\lambda} = X^{\mathsf{T}}\boldsymbol{y}$ ですから、定理 3.3 よりその正規方程式は $(X^{\mathsf{T}}X+\lambda I_p)\boldsymbol{w} = X^{\mathsf{T}}\boldsymbol{y}$ です。(1) と 系 5.2 より解はただ 1 つで、主張の式になります。$\lambda > 0$ である限り、例 5.3 のようにランクが落ちたデータでも答えが一意に定まる点が、正則化の効用です。

演習 8.4難

$P\in\mathbb{R}^{n\times n}$ が対称かつ冪等（$P^{\mathsf{T}} = P$、$P^2 = P$）であるとする。(1) $P$ の固有値は $0$ と $1$ に限ることを示してください。(2) $\operatorname{tr}P = \operatorname{rank}P$ を示してください。(3) 命題 4.3 のハット行列に対して、これが $\operatorname{tr}P = p$ と整合することを確かめてください。

解答

(1)。 $P\boldsymbol{v} = \mu\boldsymbol{v}$、$\boldsymbol{v}\ne\boldsymbol{0}$ とします。両辺に $P$ を掛けると $P^2\boldsymbol{v} = \mu P\boldsymbol{v} = \mu^2\boldsymbol{v}$ ですが、$P^2 = P$ より左辺は $P\boldsymbol{v} = \mu\boldsymbol{v}$ です。よって $(\mu^2-\mu)\boldsymbol{v} = \boldsymbol{0}$、$\boldsymbol{v}\ne\boldsymbol{0}$ より $\mu^2 = \mu$、すなわち $\mu\in\{0,1\}$ です。

(2)。 $P$ は実対称なのでスペクトル定理により直交行列 $U$ で $P = U\Lambda U^{\mathsf{T}}$、$\Lambda = \operatorname{diag}(\mu_1,\dots,\mu_n)$ と対角化できます。(1) より各 $\mu_i$ は $0$ か $1$ です。$\operatorname{tr}$ の巡回性から $\operatorname{tr}P = \operatorname{tr}(\Lambda U^{\mathsf{T}}U) = \operatorname{tr}\Lambda = \#\{i : \mu_i = 1\}$ です。一方 $U$ は正則なので $\operatorname{rank}P = \operatorname{rank}\Lambda$ で、これも $1$ の個数です。よって両者は等しくなります。

(3)。 命題 4.3 の $P = X(X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}$ は同命題の 1 より対称冪等で、その像は $\operatorname{Im}X$（同 2）ですから $\operatorname{rank}P = \dim\operatorname{Im}X = \operatorname{rank}X = p$ です。(2) と合わせて $\operatorname{tr}P = p$ となり、同命題の 4 で巡回性から直接示した結果と一致します。例 4.5 では $n = 3$、$p = 2$ で、固有値は $1, 1, 0$、トレースは $2$ でした。

参考文献

教科書（線形代数と幾何）. G. Strang, Introduction to Linear Algebra, 5th ed., Wellesley–Cambridge Press, 2016 — 第 4 章（直交性・射影・最小二乗法）。射影行列と正規方程式を幾何から説明する定番です。金谷健一『これなら分かる最適化数学：基礎原理から計算手法まで』共立出版、2005 — 最小二乗法と特異値分解を扱う章。日本語で読める入門としてまとまっています。

数値計算. G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed., Johns Hopkins University Press, 2013 — 第 5 章（直交化と最小二乗問題）。QR 分解による解法と条件数の議論はここが標準的な出典です。Å. Björck, Numerical Methods for Least Squares Problems, SIAM, 1996 — 最小二乗問題の数値解法を網羅した専門書。

統計・機械学習への接続. T. Hastie, R. Tibshirani, and J. Friedman, The Elements of Statistical Learning, 2nd ed., Springer, 2009 — 第 3 章（回帰の線形手法）。リッジ回帰・多重共線性・変数選択が本記事の枠組みの延長で説明されています。著者による PDF が公開されています。C. M. Bishop, Pattern Recognition and Machine Learning, Springer, 2006 — 第 3 章（回帰の線形モデル）。最小二乗法を最尤推定として読み直す視点が丁寧です。

歴史. S. M. Stigler, The History of Statistics: The Measurement of Uncertainty before 1900, Harvard University Press, 1986 — 第 1 部。ルジャンドルとガウスの優先権論争と、最小二乗法が受け入れられていく過程が資料に基づいて丁寧に追跡されています。

Appendix: 擬似逆行列と最小ノルム解

ランクが落ちても答えを 1 つに決める. 例 5.3 で見たように、$\operatorname{rank}X < p$ のとき最小二乗解は無数にあります。数値計算ライブラリはそのなかから 1 つを返さなければなりません。標準的な選び方は「ノルムが最小のもの」です。

命題 8.5（最小ノルム最小二乗解）✓

$X\in\mathbb{R}^{n\times p}$、$\boldsymbol{y}\in\mathbb{R}^n$ とし、最小二乗解の集合を $S$ とする。このとき $S \cap \operatorname{Im}(X^{\mathsf{T}})$ はただ 1 点からなり、その元 $\hat{\boldsymbol{w}}^{+}$ は $S$ のなかでノルムを最小にする唯一の元である。

証明（命題 8.5）

定理 4.2 より $S\ne\varnothing$ で、定理 3.3 より $S = \boldsymbol{w}_0+\operatorname{Ker}X$（$\boldsymbol{w}_0\in S$ は任意の 1 つ）です。補題 4.1 を $X^{\mathsf{T}}$ に適用すると $(\operatorname{Im}X^{\mathsf{T}})^{\perp} = \operatorname{Ker}X$ となるので、直交分解

$$\mathbb{R}^p = \operatorname{Im}(X^{\mathsf{T}}) \oplus \operatorname{Ker}X$$

が成り立ちます。これに従って $\boldsymbol{w}_0 = \boldsymbol{u}+\boldsymbol{z}$（$\boldsymbol{u}\in\operatorname{Im}X^{\mathsf{T}}$、$\boldsymbol{z}\in\operatorname{Ker}X$）と分解すると、$\boldsymbol{u} = \boldsymbol{w}_0-\boldsymbol{z} \in S$ かつ $\boldsymbol{u}\in\operatorname{Im}X^{\mathsf{T}}$ です。したがって $S\cap\operatorname{Im}(X^{\mathsf{T}})\ne\varnothing$ です。2 つの元 $\boldsymbol{u},\boldsymbol{u}'$ がこの共通部分に属せば、$\boldsymbol{u}-\boldsymbol{u}'\in\operatorname{Ker}X$（$S$ の形から）かつ $\boldsymbol{u}-\boldsymbol{u}'\in\operatorname{Im}X^{\mathsf{T}}$ ですが、直交分解の 2 つの成分の共通部分は $\{\boldsymbol{0}\}$ なので $\boldsymbol{u} = \boldsymbol{u}'$、すなわち 1 点です。これを $\hat{\boldsymbol{w}}^{+}$ と書きます。

最小性を示します。任意の $\boldsymbol{w}\in S$ は $\boldsymbol{w} = \hat{\boldsymbol{w}}^{+}+\boldsymbol{v}$（$\boldsymbol{v}\in\operatorname{Ker}X$）と書け、$\hat{\boldsymbol{w}}^{+}\perp\boldsymbol{v}$ ですからピタゴラスの定理より

$$\lVert\boldsymbol{w}\rVert^2 = \lVert\hat{\boldsymbol{w}}^{+}\rVert^2 + \lVert\boldsymbol{v}\rVert^2 \ge \lVert\hat{\boldsymbol{w}}^{+}\rVert^2$$

であり、等号は $\boldsymbol{v} = \boldsymbol{0}$ すなわち $\boldsymbol{w} = \hat{\boldsymbol{w}}^{+}$ のときに限ります。

∎

擬似逆行列. この $\hat{\boldsymbol{w}}^{+}$ は、$X$ の特異値分解 $X = U\Sigma V^{\mathsf{T}}$ から作られるムーア・ペンローズ擬似逆行列を使って書けます。

定義 8.6（ムーア・ペンローズ擬似逆行列）↗

$X\in\mathbb{R}^{n\times p}$ の特異値分解を $X = U\Sigma V^{\mathsf{T}}$（$U, V$ は直交行列、$\Sigma$ の対角成分が特異値 $\sigma_1\ge\cdots\ge\sigma_r > 0$、それ以外は $0$）とする。$\Sigma^{+}$ を、$\Sigma$ の $0$ でない対角成分を逆数に置き換えて転置した行列とし、

$$X^{+} = V\Sigma^{+}U^{\mathsf{T}} \in \mathbb{R}^{p\times n}$$

を $X$ の擬似逆行列という。$\operatorname{rank}X = p$ のときは $X^{+} = (X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}$ と一致する。

$\hat{\boldsymbol{w}}^{+} = X^{+}\boldsymbol{y}$ が成り立ちます（$X^{+}\boldsymbol{y}$ が正規方程式を満たすことと $X^{+}\boldsymbol{y}\in\operatorname{Im}X^{\mathsf{T}}$ であることを $U,\Sigma,V$ で書けば確認できます。詳しくは Golub–Van Loan の第 5 章を参照してください）。NumPy の np.linalg.lstsq や np.linalg.pinv が返すのはこの解です。例 5.3 の解集合 $(0.5, 1.5-2t, t)$ でノルムの二乗

$$0.25 + (1.5-2t)^2 + t^2 = 5t^2 - 6t + 2.5$$

を最小にするのは $t = 3/5 = 0.6$ で、$\hat{\boldsymbol{w}}^{+} = (0.5,\ 0.3,\ 0.6)^{\mathsf{T}}$ です。実際この解は核の生成元 $(0,2,-1)^{\mathsf{T}}$ と直交します（$2\cdot 0.3 - 0.6 = 0$）。ノルムの二乗は $0.7$ で、$t = 0$ の解 $(0.5, 1.5, 0)$ の $2.5$ より小さくなっています。