3.2 线性回归模型和最小二乘

如第二章所介绍，我们有一个输入变量向量 $X^T=(X_1,X_2,\dots,X_p)$，目的是预测一个实数取值的输出变量 $Y$。线性模型的表达式为：

线性模型假设回归函数 $\operatorname{E}(Y|X)$ 是线性的，或者至少足够近似于线性的。式中 $\beta_j$ 是未知的参数或系数，输入变量 $X_j$ 可以有多种类型：

• 数值类型的输入变量；
• 数值输入变量的某种变换，比如取对数、平方根、平方；
• 构建多项式表达式的基函数展开，比如 $X_2=X_1^2$，$X_3=X_1^3$。
• 分类输入变量的哑变量（dummy）或数值编码。例如，对一个有五种类别的分类变量 $G$，可构建五个变量 $X_j,j=1,\dots,5$，满足 $X_j=I(G=j)$。在 $\sum_{j=1}^5X_j\beta_j$ 中只能有一个 $X_j$ 为 1，其他都为 0，因此这组变量 $X_j$ 通过一组分类别常数就可表达 $G$ 变量的影响。
• 输入变量之间的交叉项，例如 $X_3=X_1\cdot X_2$。

不管 $X_j$ 的形式如何，模型对参数（的形式）是线性的。

我们通常有一组用来估计参数 $\beta$ 的训练数据 $(x_1,y_1),\dots,(x_N,y_N)$。每个 $x_i=(x_{i1},x_{i2},\dots,x_{ip})^T$ 是样本 $i$ 的输入向量的取值。最常见的估计方法是 最小二乘（least squares） 法，通过对残差平方和（residual sum of squares）最小化来确定系数 $\beta=(\beta_0,\beta_1,\dots,\beta_p)^T$：

从统计学的角度看，如果训练样本 $(x_i,y_i)$ 是从样本空间中独立的随机采样，那么这个目标函数的选择是合理的。既使 $x_i$ 并不是完全随机采样的，只要输出变量 $y_i$ 在条件于输入变量 $x_i$ 下相互独立，这仍是个合理的目标函数。

图 3.1 在训练样本点 $(X,Y)$ 所在的 $\mathbb{R}^{p+1}$ 空间上，演示了最小二乘法拟合的几何含义。注意式 3.2 没有任何关于模型 3.1 正确与否的假设，只是寻找训练样本数据的最佳线性拟合。不管真实的数据是怎样的，最小二乘的拟合从直觉上也是合理的；目标函数（式 3.2）衡量了（拟合线与样本的）平均欠拟合程度。

那么如何对式 3.2 最小化？记 $\mathbf{X}$ 为代表训练样本的 $N\times(p+1)$ 的矩阵，其每行对应着一个输入向量（向量第一个元素为常数 1），记 $\mathbf{y}$ 为代表训练样本输出变量的长度为 $N$ 的向量。则残差平方和可写为：

这是 $p+1$ 个参数的二次函数。对 $\beta$ 求导：

先（暂时）假设 $\mathbf{X}$ 为列满秩矩阵，则 $\mathbf{X}^T\mathbf{X}$ 是正定矩阵¹（positive definite），将一阶导数设置为 0：

可得到唯一解：

在一个变量向量 $x_0$ 点的预测值则为 $\hat{f}(x_0)=(1,x_0)^T\hat{\beta}$；训练样本的拟合值为：

或者 $\hat{y_i}=\hat{f}(x_i)$。式 3.7 中，$\mathbf{y}$ 与矩阵 $\mathbf{H}=\mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T$ 相乘后就戴上了“帽子”（$\hat{\mathbf{y}}$），所以这个矩阵也称为“帽子”矩阵。

图 3.2 在 $\mathbb{R}^N$ 空间上用另一个方式展示最小二乘估计的几何表现。将 $\mathbf{X}$ 的列向量记为 $\mathbf{x}_0,\mathbf{x}_1,\dots,\mathbf{x}_p$，其中 $\mathbf{x}_0\equiv 1$。在下面的推导中，第一个常数向量的处理方式会与其他变量相同。这些向量会生成一个 $\mathbb{R}^N$ 上的子空间，也被称为 $\mathbf{X}$ 的列空间。对 $\operatorname{RSS}(\beta)=|\mathbf{y}-\mathbf{X}\beta|^2$ 的最小化过程，等价于寻找使残差向量 $\mathbf{y}-\hat{\mathbf{y}}$ 与这个列空间正交²的 $\hat{\beta}$。式 3.5 体现了这个正交的关系，因此拟合的结果 $\hat{\mathbf{y}}$ 是 $\mathbf{y}$ 在这个列空间上的 正交投影（orthogonal projection）。帽子矩阵 $\mathbf{H}$ 实现了这个正交投影的转换计算，因此也被称为 投影矩阵。

$\mathbf{X}$ 的列向量有可能不是线性独立，那么 $\mathbf{X}$ 便不是列满秩矩阵。两个完全相关的输入变量（比如 $\mathbf{x}_2=3\mathbf{x}_1$）就是一个例子。这时 $\mathbf{X}^T\mathbf{X}$ 是个奇异（singular）矩阵，最小二乘无法得到唯一解 $\hat{\beta}$。然而拟合值 $\hat{\mathbf{y}}=\mathbf{X}\hat{\beta}$ 仍然是 $\mathbf{y}$ 在 $\mathbf{X}$ 列空间上的正交投影。只是这时有不止一种可表达这个投影的 $\mathbf{X}$ 列向量的组合方式。出现非满秩输入变量矩阵的一个最常见原因是一个或多个分类变量的编码方式出现了冗余。因此对无唯一解的自然解决方案是重新编码或者去掉 $\mathbf{X}$ 中冗余的列。大多数回归软件程序会检测 $\mathbf{X}$ 中的冗余列，并自动地采用某些去除这些冗余的策略。另一个可造成非满秩的原因是输入变量个数 $p$ 大于训练样本个数 $N$，这在信号和图像识别领域常出现。通常解决方法是过滤掉一些输入变量，或者用正则化来控制拟合（见第 5.2.3 节和第十八章）。

上述讨论没有对数据的真实概率分布做太多假设。然而为了进一推导 $\hat{\beta}$ 的性质，现假设输出变量观测样本 $y_i$ 互不相关，有固定的方差 $\sigma^2$，并且将 $x_i$ 视为非随机的常量。从式 3.6 可推导出最小二乘参数估计量的协方差矩阵为：

通常用下式来估计方差 $\sigma^2$：

注意其中的除数为 $N-p-1$ 而不是 $N$，这使得估计量 $\hat{\sigma}^2$ 是无偏估计量：$\operatorname{E}(\hat{\sigma}^2)=\sigma^2$。

为了对模型和参数做统计推断，需要做额外的假设。首先，假设式 3.1 为输出变量期望的正确模型，即 $Y$ 的条件期望对 $X_1$，……，$X_p$ 呈线性关系。其次，假设 $Y$ 与其期望的偏离是加性的（additive）且服从高斯分布（Gaussian）。因此：

其中误差项 $\varepsilon$ 为高斯随机变量，期望为 0 方差为 $\sigma^2$，即 $\varepsilon\sim\mathcal{N}(0,\sigma^2)$。

从式 3.9 可推导出：

这是一个多元正态分布，期望向量和协方差矩阵如上式所示。并且：

即自由度为 $N-p-1$ 的卡方分布。此外，$\hat{\beta}$ 与 $\hat{\sigma}^2$ 统计意义上相互独立。根据以上的概率分布，可对参数 $\beta_j$ 进行假设检验和建立置信区间。

我们构造 标准（回归）系数（standardized (regression) coefficients） 或 Z 分数（Z-score） 来检验某个系数为零的假设 $\beta_j=0$：

式中的 $v_j$ 为矩阵 $(\mathbf{X}^T\mathbf{X})^{-1}$ 的第 j 个对角线元素。在原假设下，即 $\beta_j=0$，统计量 $z_j$ 的概率分布为 $t_{N-p-1}$（自由度为 $N-p-1$ 的 $t$ 分布）。因此，如果计算得出的 $z_j$ 绝对值过大，则会拒绝原假设。若用已知真实的 $\sigma^2$ 替换等式中的估计量 $\hat{\sigma}^2$，则 $z_j$ 服从标准正态分布。当样本量足够大时，t 分布与标准正态分布的尾部百分位非常接近，所以在实践中一般会直接使用正态分布（见图 3.3）。

我们通常会需要检验一组系数是否同时显著。例如，一个有 $k$ 个类别取值的分类变量会在模型中引入多个哑变量，为了验证这个分类变量是否可从模型中去除，需要检验所有哑变量的系数同时为 0 的假设。这可使用 $F$ 统计量：

计算这个统计量需要进行两次最小二乘拟合。首先使用所有的输入变量进行最小二乘拟合，称之为模型 1，有 $p_1+1$ 个参数；然后排除掉其中一部分输入变量，或也可视为将对应的参数值固定为 0，再次进行最小二乘拟合，称之为模型 0，参数个数为 $p_0+1$。将原模型中的 $p_1-p_0$ 个输入变量的参数固定为 0，即从模型 1 得到了 模型 0，模型 0 是模型 1 的嵌套模型。上式中的脚标 1 即为模型 1，脚标 0 即为模型 0，RSS 为相应最小二乘拟合的残差平方和。$F$ 统计量衡量的是被 $\sigma^2$ 估计量标准化之后，每多增加一个输入参数可减少的拟合残差平方和。若延续之前的高斯分布假设，并且参数数量少的模型 0 为真实模型的原假设成立，则 $F$ 统计量服从概率分布 $F_{p_1-p_0,N-p_1-1}$。可证明（练习 3.1）当原假设只涉及一个参数时，等式 3.12 中的 $z_j$ 与 F 统计量的检验结果是完全等价的。当样本量 $N$ 足够大时，概率分布 $F_{p_1-p_0,N-p_1-1}$ 的分位数逼近于概率分布 $\chi^2_{p_1-p_0}$ 的分位数。

通过等式 3.10 可推导出 $\beta_j$ 的 $1-2\alpha$ 置信区间：

其中 $z^{(1-\alpha)}$ 为标准正态分布的 $1-\alpha$ 分位数：

所以标准化的实践中默认输出 $\hat{\beta}\pm2\cdot\operatorname{se}(\hat{\beta})$ 就是源自大约 95% 水平的置信区间。即使高斯分布的假设不成立，在样本量足够大的情况下 $N\rightarrow\infty$，这个区间的范围也大致是正确的，即 $1-2\alpha$ 的置信水平。

用类似的方法可得到一组输入变量参数向量 $\beta$ 在大样本下的近似置信集合：

其中 ${\chi^2_\ell}^{(1-\alpha)}$ 为 $\ell$ 个自由度的卡方分布的 $1-\alpha$ 分位数，例如 ${\chi^2_5}^{(1-0.05)}=11.1$，${\chi^2_5}^{(1-0.1)}=9.2$。而通过这个 $\beta$ 的置信集合可以得到真实函数 $f(x)=x^T\beta$ 相应的置信区间 $\{x^T\beta|\beta\in C_\beta\}$。（练习 3.2；函数的置信区域的例子另见第 5.2.2 节的图 5.4。）

3.2.1 实例：前列腺癌症

本例子中的数据来自于 Stamey et al. (1989) 的研究。其中考察了一些即将根治前列腺切除的男性病人的前列腺特异抗原（prostate-specific antigen，PSA）水平以及一些临床指标之间的相关性。特征变量有：

• 癌肿瘤体积的对数（log cancer volume，lcavol）
• 前列腺重量的对数（log prostate weight，lweight）
• 年龄（age）
• 良性前列腺增生数量的对数（log of benign prostatic hyperplasia，lbph）
• 是否存在精囊侵袭（seminal vesicle invasion，svi）
• 囊穿透的对数（log of capsular penetration，lcp）
• 格里森评分（gleason）
• 格里森评分 4 和 5 所占百分比（pgg45）

表 3.1：前列腺癌症数据集中的自变量之间的相关性矩阵。

表 3.1 为这些自变量的相关系数矩阵，可见它们之间存在着一些强相关性。第一章中图 1.1（书中第 3 页）为这些变量的两两散点图矩阵。其中 svi 为二元分类变量，gleason 为有序分类变量。从图中可看出变量 lcavol 和 lcp 都与输出变量 lpsa 有强相关性，同时两者之间也存在强相关性。通过利用这些自变量的最小二乘拟合，可以独立出每个自变量对因变量的影响。

我们将自变量的方差标准化为 1 后，对输出变量 lpsa（PSA 水平的对数） 拟合线性模型。数据集被随机分为两部分，训练集包含 67 个样本，测试集包含 30 个样本。表 3.2 为最小二乘法拟合出的参数估计值、标准误差和 Z 分数。

表 3.2：前列腺数据集的线性模型拟合结果。Z 分数按照式 3.12 计算，即参数的估计值除以其标准误差。大致上绝对值大于 2 的 Z 分数代表着这个参数在 $p=0.05$ 的置信水平下显著不为 0。

Z 分数如式 3.12 中所定义，代表着对应输入变量的显著性水平，或可理解为去掉这个变量对拟合的影响大小。大于 2 的 Z 分数即代表着大致 5% 置信水平下的显著性。（在这个数据集中，拟合中有 9 个参数，$t_{67-9}$ 分布的 1-0.025 分位数为 2.002，非常接近于 2。）自变量中 lcavol 的显著性最强，另外 lweight 和 svi 也有比较大的 Z 分数绝对值。值得注意，当 lcavol 与 lcp 同时在模型中时，lcp 并不显著。（但模型中若只有 lcp 而没有 lcavol，则 lcp 有强显著性。）

用 $F$ 统计量（式 3.13）来检验同时排除一组输入变量的假设。假如要检验是否可以排除表 3.2 中所有不显著的变量，即 age，lcp，gleason，和 pgg45，相应的 $F$ 统计量为：

其对应的 p 值为 0.17（$\operatorname{Pr}(F_{4,58}>1.67) = 0.17$），因此不显著。结论是不能拒绝这几个变量的参数同时为 0 的假设。

在测试集上的平均预测误差为 0.521。在用训练集的输出变量 lpsa 的均值作为测试集的预测时，相应的平均预测误差为 1.057，这也被称为“基础误差率（base error rate）”。因此，线性模型的误差率比基础误差率降低了超 50%。稍后章节我们会回到这个实例来对比不同的特征选择和 收缩（shrinkage） 方法。

3.2.2 高斯－马尔可夫定理

在参数 $\beta$ 的所有线性无偏估计中，最小二乘估计的方差是最小的，这是统计学中最著名的结论之一。我们将明确地表述这个结论，同时也会指出将估计量限制在无偏的条件内可能不是最优的。这引出了在本章后续会介绍的有偏估计，比如岭回归。

现在考虑对参数任意线性组合 $\theta=a^T\beta$ 的估计；例如线性模型的预测 $f(x_0)=x_0^T\beta$ 即为这样的一个线性组合。那么 $a^T\beta$ 的最小二乘估计量为：

假设 $\mathbf{X}$ 可视为常数，那么上式为输出变量向量 $\mathbf{y}$ 的一个线性函数 $\mathbf{c}^T_0\mathbf{y}$。若线性模型的假设是正确的，那么 $a^T\hat{\beta}$ 为无偏估计：

高斯-马尔可夫定理（Gauss-Markov theorem） 的内容是，其他任意一个对 $a^T\beta$ 的无偏线性估计量 $\tilde{\theta}=\mathbf{c}^T\mathbf{y}$，即 $\operatorname{E}(\mathbf{c}^T\mathbf{y})=a^T\beta$，都会满足：

上式的证明（练习 3.3）利用了三角不等式。简单起见，上述过程只考虑了对单一数值 $a^T\beta$ 的估计，在一些附加的定义后，同样的结论可扩展到针对参数向量 $\beta$（练习 3.3）。

下面考虑 $\theta$ 的估计 $\tilde{\theta}$ 的均方误差：

式中第一项为方差，第二项为平方偏差。高斯-马尔可夫定理说明在所有的线性无偏估计量中，最小二乘估计量的均方误差是最小的。然而根据上式，如果能够以较小偏差的代价带来较大的方差降低，则有可能存在一个均方误差更小的有偏的估计。有偏估计是很常用的。所有将最小二乘法中的部分参数向零缩减（shrink）或直接设为零的方法都可产生有偏估计。本章稍后介绍的变量选择和岭回归方法都是有偏估计的例子。在实际应用场景中，模型都是对现实的简化，因此得到的都是有偏估计；模型选择就是在偏差和方差之间寻找合适的平衡点。第七章会更详细讨论相关内容。

在第二章提到过，均方误差是与预测的准确程度紧密关联的指标。假设要对在输入变量 $x_0$ 处的值做预测：

$$Y_0 = f(x_0) + \varepsilon_0 \tag{3.21}$$

则某个估计量 $\tilde{f}(x_0) = x_0^T\beta$ 的期望预测误差为

$$\begin{align} E(Y_0 - \tilde{f}(x_0))^2 &= \sigma^2 + E(x_0^T\tilde{\beta} - f(x_0))^2 \\ &= \sigma^2 + \operatorname{MSE}(\tilde{f}(x_0)) \tag{3.22} \end{align}$$

因此期望预测误差与均方误差仅相差一个代表观测本身的随机性的方差常数 $\sigma^2$。

3.2.3 从简单一元回归到多元回归

当 $p>1$ 时，式 3.1 中的线性模型称为多元线性回归模型。本节从单变量线性模型（$p=1$）入手，以此能够更好地理解多元的最小二乘解（式 3.6）。

先假设一个无常数项的单变量模型：

则最小二乘的估计和残差为：

定义 $\mathbf{y}=(y_1,\dots,y_N)^T$，$\mathbf{x}=(x_1,\dots,x_N)^T$，以及两个向量的内积³为：

则式 3.24 可用向量和内积来表达：

这个简单的单变量回归是多元线性回归的基础。下面假设多元线性回归的多个输入变量向量 $\mathbf{x}_1$，$\mathbf{x}_2$，……，$\mathbf{x}_p$（矩阵 $\mathbf{X}$ 的列向量）彼此正交；即对于任意的 $j\neq k$，有 $\langle\mathbf{x}_j,\mathbf{x}_k\rangle=0$。则不难推导出多元最小二乘的估计 $\hat{\beta}_j$ 为 $\langle\mathbf{x}_j,\mathbf{y}\rangle/\langle\mathbf{x}_j,\mathbf{x}_j\rangle$，即这个变量的单变量回归估计。也就是说当输入变量彼此正交时，它们在模型参数的估计中彼此没有影响⁴。

通常在精心设计的实验中可以得到正交的输入变量，但在观测数据集中几乎不可能存在。因此我们需要通过对输入变量的正交化才可使用上述的结果。接下来考虑一个常数项和一个变量 $\mathbf{x}$ 作为输入变量的模型。则 $\mathbf{x}$ 的最小二乘系数估计为：

其中的 $\bar{x}=\sum_ix_i/N$，$\mathbf{1}=\mathbf{x}_0$ 为一个长度为 $N$ 每个元素为 1 的向量。式 3.27 中的估计可被分解为使用了两次简单回归（式 3.26）：

1. $\mathbf{x}$ 对常数向量 $\mathbf{1}$ 做回归，生成残差 $\mathbf{z}=\mathbf{x}-\bar{x}\mathbf{1}$。
2. $\mathbf{y}$ 对残差向量 $\mathbf{z}$ 做回归，得到系数的估计 $\hat{\beta_1}$。

其中的“$\mathbf{b}$ 对 $\mathbf{a}$ 做回归（regress $\mathbf{b}$ on $\mathbf{a}$）”指的是以 $\mathbf{b}$ 为因变量 $\mathbf{a}$ 为自变量的无截距项的单变量回归，其斜率的估计为 $\hat{\gamma}=\langle\mathbf{a},\mathbf{b}\rangle/\langle\mathbf{a},\mathbf{a}\rangle$，残差向量为 $\mathbf{b}-\hat{\gamma}\mathbf{a}$。这个过程称之为用 $\mathbf{a}$ 对 $\mathbf{b}$ 进行调整、或正交化（orthogonalize）。

第一步中用常数向量 $\mathbf{x}_0=\mathbf{1}$ 正交化 $\mathbf{x}$。第二步中用第一步的残差 $\mathbf{z}$ 为输入变量进行简单单变量回归⁵。图 3.4 以两个输入变量为例演示了这个过程，正交化不会改变输入向量 $\mathbf{x}_1$ 和 $\mathbf{x}_2$ 生成的子空间，而是生成了这个子空间的一组正交基。

算法 3.1 将上述流程推广至 $p$ 个输入变量的回归中。注意其中第二步中生成的残差向量$\mathbf{z}_0$，……，$\mathbf{z}_{j-1}$ 彼此正交，所以其中计算出的单变量回归的系数（$\hat{\gamma_{lj}}$）实际上也是多元回归的系数。

算法 3.1：逐步正交化（successive orthogonalization）回归

1. 初始化常数向量 $\mathbf{z}_0=\mathbf{x}_0=\mathbf{1}$.
2. 对 $j=1,2,\dots,p$ 依次进行：
   通过 $\mathbf{x}_j$ 对 $\mathbf{z}_0$，$\mathbf{z}_1$，……，$\mathbf{z}_{j-1}$ 的回归，得到系数估计 $$\hat{\gamma}_{\ell j} = \frac {\langle \mathbf{z}_\ell, \mathbf{x}_j \rangle} {\langle \mathbf{z}_\ell, \mathbf{z}_\ell \rangle}, \ell = 0, \dots, j−1$$ 以及回归的残差 $$\mathbf{z}_j = \mathbf{x}_j - \sum_{k=0}^{j-1} \hat{\gamma}_{kj}\mathbf{z}_k$$
3. 通过 $\mathbf{y}$ 对残差 $\mathbf{z}_p$ 的回归得到参数估计 $\hat{\beta}_p$。

算法最终的结果为多元回归中的第 $p$ 个参数的估计：

算法第二步中的逐步正交过程，使得每个向量 $\mathbf{x}_j$ 均为一组向量 $\mathbf{z}_k$（$k\leq j$）的线性组合。所有的 $\mathbf{z}_j$ 彼此正交，因此它们构成了 $\mathbf{X}$ 的列空间的一组基（basis）。最小二乘拟合 $\hat{\mathbf{y}}$ 是输出变量向量在这个空间上的投影，所以它也是这组基的线性组合。由于基中只有 $\mathbf{z}_p$ 中包含 $\mathbf{x}_p$（并且系数为 1），其在投影的正交基线性组合中系数为式 3.28，因此这个系数也是原多元回归中 $\mathbf{x}_p$ 的最小二乘系数。通过改变输入变量的脚标顺序，将不同的输入变量放在最后一位，则重复算法 3.1 就可以得到所有参数估计。

算法 3.1 不仅提供了另一个计算最小二乘的方法，同时也揭示了多元回归中存在相关性的输入变量之间的影响。这个过程说明了多元回归中的第 j 个回归系数，等价于 $\mathbf{y}$ 对 $\mathbf{x}_{j\cdot012\dots(j-1)(j+1)\dots p}$（$\mathbf{x}_j$ 对所有其他输入变量 $\mathbf{x}_0$，$\mathbf{x}_1$，……，$\mathbf{x}_{j-1}$，$\mathbf{x}_{j+1}$，……，$\mathbf{x}_p$ 回归的残差）的单变量回归系数。

直观的来说：多元回归系数 $\hat{\beta}_j$ 体现的是排除了其他变量的影响后，$\mathbf{x}_j$ 对 $\mathbf{y}$ 的影响。

若 $\mathbf{x}_p$ 与其他一个或几个输入变量有强相关性，那么最终的残差向量 $\mathbf{z}_p$ 会十分接近于 0，通过式 3.28 可知其系数估计 $\hat{\beta}_p$ 不稳定。与之强相关的其他输入变量也同样如此。此时可能出现的是，这组变量的 Z 评分都不大（如图 3.2），其中每个变量都可以被排除，然而我们却无法将这组变量全部排除。从式 3.28 也可得到方差估计（式 3.8）的另一个表达式：

也就是说，对 $\hat{\beta}_p$ 估计的准确度取决于残差向量 $\mathbf{z}_p$ 的长度，即在 $\mathbf{x}_p$ 的信息中有多少是无法被其他输入变量解释的。

算法 3.1 也被称为多元回归的 格拉姆-施密特正交化（Gram–Schmidt procedure），它也是一个实践中最小二乘估计的数值计算方法。算法最终可获得参数估计 $\hat{\beta}_p$，以及最小二乘的拟合结果（练习 3.4）。

算法 3.1 中的第二步可以用矩阵的形式表达：

其中 $\mathbf{Z}$ 的列（依次）为 $\mathbf{z}_j$，$\mathbf{\Gamma}$ 为上三角矩阵，矩阵中的元素为 $\hat{\gamma}_{kj}$。定义对角矩阵 $\mathbf{D}$，其第 j 个对角线元素为 $\mathbf{D}_{jj}=|\mathbf{z}_j|$，则：

便是所谓的矩阵 $\mathbf{X}$ 的 QR 分解（QR decomposition）。其中的 $\mathbf{Q}$ 为 $N\times(p+1)$ 的正交矩阵，且 $\mathbf{Q}^T\mathbf{Q}=\mathbf{I}$，$\mathbf{R}$ 为 $(p+1)\times(p+1)$ 的上三角矩阵。QR 分解为矩阵 $\mathbf{X}$ 的列向量空间提供了一个方便的正交基。例如，最小二乘的解可以很方便的写为：

由于 $\mathbf{R}$ 为上三角矩阵，其逆矩阵很容易获得（练习 3.4）。

3.2.4 多输出变量回归

假设有多个输出变量 $Y_1$，$Y_2$，……，$Y_K$，都用 $p$ 个输入变量 $X_0$，$X_1$，$X_2$，……，$X_p$ 来预测。则对每一个输出变量，都有线性模型：

假设样本大小为 $N$，可以将模型写成矩阵形式：

$$\mathbf{Y} = \mathbf{X}\mathbf{B} + \mathbf{E} \tag{3.36}$$

其中 $\mathbf{Y}$ 为一个 $N\times K$ 的因变量矩阵，其中第 $ik$ 位置的元素为 $y_{ik}$；$\mathbf{X}$ 为 $N \times(p+1)$ 的输入变量矩阵；$\mathbf{B}$ 为 $(p+1)\times K$ 的参数矩阵；$\mathbf{E}$ 为 $N\times K$ 的误差项矩阵。对单输出变量的损失函数（式 3.2）的简单推广可得到：

最小二乘法的估计量和单输出变量时的解一致，但其中的矩阵维度不同：

因此第 k 个输出变量的系数估计等同于 $\mathbf{x}_0$，$\mathbf{x}_1$，……，$\mathbf{x}_p$ 对 $\mathbf{y}_k$ 的最小二乘回归。多个输出变量不会影响彼此的最小二乘估计。

但如果式 3.34 中的误差项 $\varepsilon=(\varepsilon_1,\dots,\varepsilon_K)$ 彼此存在相关，则可能需要将式 3.37 的定义针对多输出变量做调整。假设 $\operatorname{Cov}(\varepsilon)=\mathbf{\Sigma}$，则多元加权判别标准为：

其与多元高斯分布的等式比较相像。$f(x)$ 的取值为一个向量 $(f_1(x),\dots,f_K(x))^T$，$y_i$ 为观测样本 i 的长度为 K 的向量。然而推导得出的解与式 3.39 一致，即误差无相关性时的分别进行的 K 个回归解（练习 3.11）。在样本间有不同的协方差矩阵 $\mathbf{\Sigma}_i$ 时，以上结论不再成立，输出变量之间会有相互影响。

第 3.7 节会继续讨论多个输出变量的问题，以及一些联合回归可改善预测结果的场景。

本节练习

练习 3.1

证明在检验是否可从模型中排除单个系数时，$F$ 统计量（式 3.13）等于相应的 $Z$ 分数（式 3.12）的平方。

练习 3.2

给定两个变量 $X$ 和 $Y$ 的数据样本，拟合一个三次多项式回归模型 $f(X)=\sum_{j=0}^3\beta_j X^j$。除拟合曲线外，你需要围绕曲线的 95% 水平置信区间带。考虑以下两种方法：

1. 在每个点 $x_0$，为线性函数 $a^T\beta=\sum_{j=0}^3\beta_jX^j$ 构建 95% 置信区间。
2. 先如式 3.15 构建 $\beta$ 的 95% 置信区间，再从中生成函数 $f(x_0)$ 的置信区间。

这两种方法有何区别？哪个更可能得出更宽的置信区间带？通过一个小模拟试验来比较两种方法。

练习 3.3

高斯-马尔可夫定理（Gauss-Markov theorem）：

1. 证明高斯－马尔可夫定理：一个参数 $a^T\beta$ 的最小二乘估计的方差不大于任意其他 $a^T\beta$ 的线性无偏估计的方差。
2. 矩阵不等式 $\mathbf{B}\preceq\mathbf{A}$ 的定义为 $\mathbf{A}-\mathbf{B}$ 是半正定矩阵。若 $\hat{\mathbf{V}}$ 为 $\beta$ 最小二乘估计的协方差矩阵，$\tilde{\mathbf{V}}$ 为任意其他线性无偏估计的协方差矩阵，则 $\hat{\mathbf{V}}\preceq\tilde{\mathbf{V}}$。

练习 3.4

利用 $\mathbf{X}$ 的 QR 分解，试说明在一轮格拉姆-施密特（Gram－Schmidt）正交化流程（算法 3.1）后即可得到最小二乘系数向量的估计。

练习 3.11

证明多输出变量线性回归问题（式 3.40）的解为式 3.39。如果每个观测样本的协方差矩阵 $\mathbf{\Sigma}_i$ 不同，会发生什么变化？