http://www.ma-xy.com
目录
第一章 回归模型 1
1.1 问题说明 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 参数回归 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2.1 线性回归 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2.2 广义线性回归 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2.3 参数估计 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2.4
正则化
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.2.5 随机梯度下降算法及其变体 . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.3 贝叶斯回归模型 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.3.1 贝叶斯统计基础 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.3.2 贝叶斯线性回归模型 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.4 RVM 稀疏向量机 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
1.4.1 RVM 的建立 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
1.4.2 最大边缘似然方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
1.4.3 EM 算法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
1.4.4 快速序列算法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
1.4.5 核方法扩展-多核方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
1.5 MATLAB 回归简介 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
1.5.1 MATLAB 回归命令 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
1.5.2 MATLAB 回归示例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
1.6 RVM 工具箱 - SB2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
1.7 非参数回归 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
1.7.1 核估计 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
1.7.2 样条拟合 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
1.7.3 正交级数回归 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
1.7.4 小波核估计 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
1
http://www.ma-xy.com
第一章 回归模型
1.1 问题说明
在统计基础章节中，我们讨论了变量的参数估计假设检验、分类变量对连续变量的方差分析
以及变量相关性检验等问题，这些仅仅是基础统计，在实际处理数据时，我们更关心的是变量之
间的关系 (x, y 之间是什么函数关系)？下面，将要介绍变量关系模型。考虑如表 (1.1) 的数据
表 1.1: 变量关系模拟数据
data X Y
x
1
x
2
··· x
p
y
1
y
2
. . . y
q
1
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
2
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
设变量 X 含有 p 个分量，X = (x
1
, x
2
, . . . , x
p
)
T
，其中，x
i
是一个变量，如果 x
i
是随机变
量，X 就是随机向量。设共有 n 个样本，则 X 的样本数据是一个 n ×p 的矩阵数据。设变量 Y
含有 q 个分量，Y = (y
1
, y
2
, . . . , y
q
)
T
，共有 n 个样本，故其样本数据为 n × q 的矩阵。现在，要
求 X 和 Y 的关系式，即 Y = f(X) 中的映射 f。
先从简单的问题开始。上述问题是在 p − q 维上考虑的 (即 f : R
p
→ R
q
)，我们将其简化为
p = q = 1 的情况，那么 y = f(x) 就是 xOy 平面上的一个曲线。我们来考虑下面两种数据 (样
本) 情况：
¬统计数据：某社区家庭年月可支配收入与消费支出数据如表 (1.2)
¬
¬
《计量经济学 (第三版)》李子奈
1
http://www.ma-xy.com
1.1 问题说明 第一章 回归模型
表 1.2: 某社区家庭年月可支配收入与消费支出数据
每月家庭可支配收入 X 800 1100 1400 1700 2000 2300 2600 2900 3200 3500
每月家庭消费支出 Y
561 638 869 1023 1254 1408 1650 1969 2090 2299
594 748 913 1100 1309 1452 1738 1991 2134 2321
627 814 924 1144 1364 1551 1749 2046 2178 2530
638 847 979 1155 1397 1595 1804 2068 2266 2629
935 1012 1210 1408 1650 1848 2101 2354 2860
968 1045 1243 1474 1672 1881 2189 2486 2871
1078 1254 1496 1683 1925 2233 2552
1122 1298 1496 1716 1969 2244 2585
1155 1331 1562 1749 2013 2299 2640
1188 1364 1573 1771 2035 2310
1210 1408 1606 1804 2101
1430 1650 1870 2112
1485 1716 1947 2200
2002
共计 2420 4950 11495 16445 19305 23870 25025 21450 21285 15510
上表的数据可以整理为下面的一般形式
x = x
1
, x
2
, . . . , x
n
x中有重复
y = y
1
, y
2
, . . . , y
n
实验数据：
x = x
1
, x
2
, . . . , x
n
x中无重复
y = y
1
, y
2
, . . . , y
n
统计数据中允许样本取相同的值，如果这时我们要求解 y = f (x)，则用 (参数) 回归等统计
方法；实验数据中往往不允许样本取相同的值 (也不是绝对的)，如果要求 y = f(x)，则使用拟
合等方法。但是有时候二者的界限并不是很清楚，或者说回归和拟合是同一个工作。前面我们说
过，无论是微分方程、函数逼近还是回归拟合，本质工作都是寻找 f 。我们用图形来展示这两种
数据，如图 (1.1) 所示
http://www.ma-xy.com 2 http://www.ma-xy.com
http://www.ma-xy.com
第一章 回归模型 1.2 参数回归
图 1.1: 统计数据和实验数据的图示
我们的问题是：y 与 x 是什么关系？即求 f 。从局部来看问题，在给定 x = x
i
时，对应的
y
i
的取值是多少？这里，记 y 的估计值为 ˆy，即 ˆy = f (x)。在给定 f 的具体形式之后，每给一
个 x
i
，都会对应一个 ˆy
i
= f(x
i
)。当然，估计量 ˆy
i
和真实数据 y
i
之间会有误差，记误差为 ε，ε
i
是 ˆy
i
和 y
i
的误差。在模型中，有时候我们并没有把模型输出写成 ˆy，但是要知道，模型的输出
就是真实数据 y 的估计 ˆy。
1.2 参数回归
1.2.1 线性回归
上面简单介绍了我们要处理的问题，下面将介绍一些求解这类回归拟合问题的方法。与后面
的非参数回归相比，参数回归是关系式由含参 θ 的 y = f(x|θ) 决定的回归模型，我们在 x
i
点对
y
i
的估计就是用 f(x
i
|θ) 进行的。而非参数回归是依据样本数据直接给出 x
i
点的 y
i
的估计值，
不需要外来参数 θ，并且往往也不需要 y = f (x) 有具体的形式。下面，先来介绍参数回归，然后
再介绍非参数回归。
参数回归的一般形式为
y = f(x|θ) + ε
其中：θ 是未知参数，f 是含参关系式，ε 是误差项，x, y 是一维变量，x, y ∈ R。f 的形式有参
数 θ 确定，求 f 也即求 θ。
如果假设 f 为线性的，即 f (x|θ) = w
0
+ w
1
x，则有
y = w
0
+ w
1
x + ε
ε
i
iid
∼ N (0, σ
2
)
y
i
|x
i
∼ N(w
0
+ w
1
x
i
, σ
2
)
这里并不对 ε 做过多的讨论，假设 ε
i
iid
∼ N (0, σ
2
)，当然也可以假设 ε 服从其它分布。有关 ε
以及模型检验、参数检验和估计量性质等问题，在后面的 Logistics 回归中会做一些简单的介绍。
http://www.ma-xy.com 3 http://www.ma-xy.com
http://www.ma-xy.com
1.2 参数回归 第一章 回归模型
如果把 x 扩展到 p 维，X = (x
1
, x
2
, . . . , x
p
)
T
，则有多元线性回归模型
y = w
0
+ w
1
x
1
+ w
2
x
2
+ ··· + w
p
x
p
+ ε
ε
i
iid
∼ N (0, σ
2
)
其中：X ∈ R
p
，y ∈ R，ε 为误差项。如果继续把 y ∈ R 扩展到 q 维，则多元线性回归模型表示为
Y = AX + ε
ε ∼ N(0, Σ)
Y |X = N (AX, Σ)
其中：A, X, Y, ε 皆是矩阵 (向量) 形式，ε 为误差项，是一个 n ×p 维矩阵。
1.2.2 广义线性回归
上面建立的线性回归模型仅能解决 x, y 是线性关系的情况，而实际数据所呈现出来的趋势
往往并不是线性的，如图 (1.2)(a)(b) 所示
(a) 线性回归 (b) 非线性回归
图 1.2: 非线性回归示意图
为了解决图 (1.2)(b) 的非线性情况，我们考虑下面两种模型
y = ax
2
+ bx + c + ε
y = w
0
+ w
1
ϕ(x) + ε
第一种模型是一个二次多项式拟合模型，它要求预先给出的数据大致是 2 次多项式形式的，
然后再对 f 中的参数 a, b, c 进行求解，这样做的缺点是：许多情况下，并不知道数据是几次的。
第二种模型是一个带映射 ϕ 的线性模型，由于映射 ϕ 的广泛性，第二种模型可以拟合很多数据，
但是此模型的难点是 ϕ 的确定。由于第二种模型形式上仍然是线性的，我们称这种模型为广义
线性模型，后面将主要讨论这种模型，并且假设总能预先找到合理的 ϕ。
http://www.ma-xy.com 4 http://www.ma-xy.com
http://www.ma-xy.com
第一章 回归模型 1.2 参数回归
1.2.3 参数估计
上面给出了线性模型和广义线性模型，模型中有参数 w = (w
1
, w
2
, . . . ) 和 σ
2
是待求的，下
面将介绍一些参数估计方法。考虑如下 (广义) 线性模型
y = w
0
+ w
1
ϕ(x) + ε
ˆy = w
0
+ w
1
ϕ(x)
ε
i
iid
∼ N (0, σ
2
)
y
i
|x
i
∼ N(w
0
+ w
1
ϕ(x
i
), σ
2
)
(1.1)
其中：x, y ∈ R，ϕ 为一映射。有样本数据 x
i
, y
i
，样本数为 n，模型中 w
0
, w
1
, σ 为待求参数。我们的
目标是求 f，而 f 由 w
0
, w
1
, σ 决定，所以目标变为 w
0
, w
1
, σ 的参数估计问题。记 θ = (w
0
, w
1
, σ)，
对于求 θ，很自然想到：如果问题是一个优化问题就好了，我们在参数空间 θ ∈ Θ 中寻找最优参
数，使得某个给定的目标最小。最终迭代得到的 θ 即为参数 θ 的估计值，记为
ˆ
θ 或者 θ
∗
。
极大似然估计 ML
在前面的模型 (1.1) 中，给出了样本概率 y
i
|x
i
∼ N(w
0
+ w
1
ϕ(x
i
), σ
2
)，还可以写为 y
i
|x
i
∼
N(ˆy
i
, σ
2
) = N(f (x
i
), σ
2
)。由极大似然估计的思想 (即“样本出现的概率最大”) 想到，优良的 θ
应该会使样本数据出现的概率最大。由于样本是独立同分布的，可以给出样本的联合概率密度
L = P {y
1
, y
2
, . . . , y
n
} =
n
i=1
p(y
i
)
由此可以求解 w，使 L 最大，有
max
w
L(w) =
n
i=1
p(y
i
)
=
n
i=1
1
√
2πσ
2
e
−
(y
i
− ˆy
i
)
2
2σ
2
=
1
(2π)
n
2
σ
n
e
−
1
2σ
2
(y −ˆy)
′
(y −ˆy)
注：已知 ϕ，给出一个 w，则有 f(x|w)，即 ˆy
i
= f(x
i
|w)，从而有 y
i
− ˆy
i
，从而有 L(w)。这里
的 ˆy = (ˆy
1
, ˆy
2
, . . . , ˆy
n
)。
要求 L(w) 的最大值，由于 L(w) 与 ln L(w) 取最大值时的 w
∗
相同，所以 max
w
L(w) 等价于
max
w
ln L(w)，即上述优化问题变为
max
w
ln L(w) = −n ln(
√
2πσ) −
1
2σ
2
(y − ˆy
′
(y − ˆy))
并且，注意到上式右边第 1 项与 w 无关，是常数，故优化问题变为
min
w
(y − ˆy)
′
(y − ˆy) = min
w
n
i=1
(y
i
− ˆy
i
)
2
http://www.ma-xy.com 5 http://www.ma-xy.com
http://www.ma-xy.com
1.2 参数回归 第一章 回归模型
由极值原理
∂ ln L
∂w
= 0
可以求 w。
最小二乘估计 OLS
上面的极大似然估计 (ML) 的目标是求 w，使得样本出现的概率最大，ML 方法是完全基于
统计思想的 (在概率下使用)。下面，我们来介绍一种应用更普遍的参数估计方法 - 最小二乘法
(OLS)。
最小二乘法的目标是求 w，使得 y 和 ˆy 的离差平方和最小，即
min
w
Q(w) =
n
i=1
(y
i
− ˆy
i
)
2
(1.2)
其中：ˆy
i
= w
0
+ w
1
ϕ(x
i
)。其实，我们还可以设定其它目标，比如：求 w 使得离差绝对值的和最
小；求 w 使得离差最大值最小等等
min
w
n
i=1
|y
i
− ˆy
i
|
min
w
max
1⩽i⩽n
|y
i
− ˆy
i
|
min
w
p
n
i=1
(y
i
− ˆy
i
)
p
总之，目标很关键也很灵活，并且，可以看到 OLS 最终的目标和 ML 的一样。求 Q(w) 的最小
值，即 Q 对 w
0
, w
1
求导，然后令导数为 0，有
∂Q
∂w
0
= 0
∂Q
∂w
1
= 0
⇒
ˆw
0
= ¯y − ˆw
1
¯
ϕ(x)
ˆw
1
=
ϕ(x
i
)y
i
ϕ(x
i
)
2
上面的模型 (1.1) 是在二维 x, y ∈ R 中讨论的，并且没有给出参数 σ 的估计 (这个参数的估
计应该很明了)。现在，我们将 x 推广到 p 维，考虑如下广义多元线性回归模型
y = w
0
+
m−1
j=1
w
j
ϕ
j
(x) + ε
ε
i
iid
∼ N (0, σ
2
)
y
i
|x
i
∼ N(ˆy
i
, σ
2
)
http://www.ma-xy.com 6 http://www.ma-xy.com
http://www.ma-xy.com
第一章 回归模型 1.2 参数回归
将上面的模型写为矩阵的形式，有
y = ϕ(x)w + ε
当然也可以写为 y = w
T
ϕ(x)+ε，只是行列调换而已。其中：m 为映射 ϕ 的大小，ϕ = (ϕ
1
, ϕ
2
, . . . , ϕ
m
)，
ϕ(x) 是一个 n × m 维矩阵
ϕ(x) =
ϕ
0
(x
1
) ϕ
1
(x
1
) . . . ϕ
m−1
(x
1
)
ϕ
0
(x
2
) ϕ
1
(x
2
) . . . ϕ
m−1
(x
2
)
.
.
.
.
.
.
.
.
.
.
.
.
ϕ
0
(x
n
) ϕ
1
(x
n
) . . . ϕ
m−1
(x
n
)
n×m
如果取特殊的 ϕ，我们可以令 m = p，或者说每一个变量 x
j
都加一个 ϕ
j
，即
y = w
0
+
p
j=1
ϕ
j
(x
j
) + ε
我们仍然采用最小二乘来求参数 w = (w
0
, w
1
, . . . , w
m−1
)，有
min
w
Q(w) =
n
i=1
ε
2
i
= ε
′
ε
= ( y −ϕ(x)w)
′
(y − ϕ(x)w)
= ||y −ϕ(x)w||
2
由极值原理，有
∂Q
∂w
= 0
⇒
∂(y − ϕ(x)w)
′
(y − ϕ(x)w)
∂w
= 0
⇒
∂(y
′
y − w
′
ϕ(x)
′
y − y
′
ϕ(x)w + w
′
ϕ(x)
′
ϕ(x)w)
∂w
= 0
⇒
∂
∂w
(y
′
y − 2y
′
ϕ(x)w + w
′
ϕ(x)
′
ϕ(x)w) = 0
即
−ϕ(x)
′
y + ϕ(x)
′
ϕ(x)w = 0
推得
w = (ϕ(x)
′
ϕ(x))
−1
ϕ(x)
′
y (1.3)
注：其本质是这样的
y = ϕ(x)w
ϕ(x)
′
y = ϕ(x)
′
ϕ(x)w
(ϕ(x)
′
ϕ(x))
−1
ϕ(x)
′
y = w
http://www.ma-xy.com 7 http://www.ma-xy.com
http://www.ma-xy.com
1.2 参数回归 第一章 回归模型
1.2.4 正则化
在上面求得的权重 w 计算公式 (1.3) 中，w 被写为
w = (ϕ(x)
′
ϕ(x))
−1
ϕ(x)
′
y
为简便，我们将其简略为
w = (x
T
x)
−1
x
T
y (1.4)
对于计算公式 (1.4)，有两个值得讨论的地方：¬ x
T
x 是否可逆，即推导时 x
T
x 能除过去吗？
在 x
T
x 可逆的情况下，(x
T
x)
−1
的计算量如何？可以想象，当 n 越大时，计算量越大。
对于一个方阵 x
T
x 而言，如果行列式 |x
T
x| 不满秩，则 x
T
x 不可逆，其特征根中至少有一
个是 0，矩阵为奇异矩阵；如果 |x
T
x| = 0 ，则 x
T
x 可逆，即 (x
T
x)
−1
可求。
由数据集给出的 x
T
x 肯定是可逆的，但是 x
T
x 的行列式 |x
T
x| 的值可能非常小，这导致 x
T
x
的一个很小的差异，就会引起 (x
T
x)
−1
很大的波动 (可以想象为
1
a
，a 很小，则 a 的一个微小的
波动
1
a
都变化很大，从 y =
1
x
的图像可以明显看出)。即当 |x
T
x| → 0 时，w 的估计具有不稳定
性。令 A = x
T
x，则
A
−1
=
1
|A|
A
∗
其中：A
∗
是 A 的伴随矩阵。对于近似奇异的情况 |A| → 0，|A| 的一个小变化会导致
1
|A|
巨大的
变化，从而 A
−1
变化巨大。对于这种近似奇异的情况，我们只要将 |A| 调大一些就好了，由于
|A| → 0 时，|A| =
n
i=1
λ
i
，不妨将 A 修正为
A := A + λI
其中：λ 为常参数，I 为单位矩阵。
由上面的修正思想，我们可以得到如下的 w 的修正计算公式
w = (x
T
x + λI)
−1
x
T
y (1.5)
其中：λ 为 turning parameter。与 (1.5) 式对应的最小二乘优化问题为岭回归 (ridge regression)
min
w
Q(w) = ||y − xw||
2
+ λ||w||
2
=
n
i=1
(y
i
− ˆy
i
)
2
+ λ
m
j=0
w
2
j
我们称 ||w||
2
项为正则项 (罚项)。上面岭回归的正则项是二次的，当然我们可以用其它的形式，
比如 lasso regression 中用的一次罚项 (1996.Tibshirami)
min
w
Q(w) = ||y − wx||
2
+ λ||w||
1
http://www.ma-xy.com 8 http://www.ma-xy.com
http://www.ma-xy.com
第一章 回归模型 1.2 参数回归
其中：||w||
1
=
m
j=1
|w
j
|。但是，lasso 并不能像岭回归和普通线性回归那样写出显式的 w 的计
算公式。更一般的，我们还可以写出正则项为 λ||w||
p
时的最小二乘优化目标
min
w
Q(w) = ||y − ˆy||
2
+ ||w||
p
我们来观察一下不同 p 时，正则项的轮廓线图，如图 (1.3) 所示
图 1.3: 正则项的轮廓线
一般情况下，最好令 p ⩾ 1，因为当 p < 1 时，集合变为非凸，从而优化问题也变为一个非
凸优化。
上面是在 x 上讨论的，我们将 x 还原为 ϕ(x)。当 p = 2 时，ϕ(x) 对应的岭回归的最小二乘
优化目标为
min
w
Q(w) =
1
2
||y − ˆy||
2
+
1
2
λ||w||
2
=
1
2
n
i=1
(y
i
− ϕ(x
i
)w)
2
+
1
2
λw
T
w (1.6)
其 w 计算公式为
w = (ϕ
T
(x)ϕ(x) + λI)
−1
ϕ
T
(x)y
当 p = 1 时，ϕ(x) 对应的 lasso regression 的最小二乘优化目标为
min
w
Q(w) =
1
2
n
i=1
(y
i
− ϕ(x
i
)w)
2
+
1
2
λ
m
j=1
|w
j
| (1.7)
如果 λ 很大，那么某些 w
j
就会变为 0，从而产生稀疏 (Spase) 模型。值得一提的是，当 ϕ
T
ϕ 接
近奇异时，可以通过奇异值分解 (SVD) 来解决这一问题，并且记 ϕ
†
= (ϕ
T
ϕ)
−1
ϕ
T
是矩阵 ϕ 的
Moore-Penrose 伪逆矩阵。关于 ridge 和 lasso 我们后面还会介绍。
http://www.ma-xy.com 9 http://www.ma-xy.com
http://www.ma-xy.com
1.2 参数回归 第一章 回归模型
1.2.5 随机梯度下降算法及其变体
上面建立的回归模型如：广义线性回归 (1.2)、岭回归 (1.6) 和 lasso 回归 (1.7) 都有相应的
优化算法。我们写出最小二乘问题一般形式
min
θ
J(θ) =
1
n
f
i
(θ) (1.8)
其中：f
i
就相当于最小二乘中的 (y
i
− ˆy
i
(w))
2
，我们称其为样本点 (x
i
, y
i
) 的损失函数。注意：这
里的一般形式 (1.8) 并没有包含正则化项。MATLAB 最优化工具箱使用两种算法来求解有约束
线性最小二乘问题：
1. 有效集算法用于求解有边界的问题以及线性不等式或等式。
2. 信赖域反射用于求解只带有边界约束的大规模问题。
该工具箱使用两种算法来求解非线性最小二乘问题：
1. 信赖域反射算法采用信赖域方法来实现 Levenberg-Marquardt 算法。它用于无约束和有边
界约束的问题。
2. Levenberg-Marquardt 算法实现标准的 Levenberg-Marquardt 算法。它用于无约束问题。
MATLAB 求解最小二乘优化问题的命令为 lsqcurvet，应用示例
1 x = Data ( : , 1 ) ; y = Data ( : , 2 ) ;
2 F = @(w, xdata )w(1) *exp(−w( 2 ) *xdata ) + w( 3 ) *exp(−w( 4) *xdata ) ;
3 w0 = [ 1 1 1 0 ] ;
4 [w, resnorm , ~ , e x i t f l a g , output ] = l s q c u r v e f i t (F, w0, x , y )
5 pl o t ( x , y , ’ ro ’ ) , hold on
6 pl o t ( x , F(w, x) ) , hold o f f
7 t i t l e ( ’ 非 线 性 最 小 二 乘 ’ )
8
关于最小二乘优化问题的一般性解法，可以参考最优化部分相应的章节。lsqcurvet 采用的方法
都是确定性算法：每一迭代步长是多少，迭代方向是哪里，都是完全确定的。下面，我们来介绍
机器学习模型中常用的优化方法 - 随机梯度下降方法及其变体，这种方法具有随机性，可以跳出
局部极小值点。
为了和一般的参考文献统一，我们将最小二乘模型中的符号变为下列形式
J(θ) =
1
2n
n
i=1
y
(i)
− h
θ
x
(i)
2
其中：y
(i)
为第 i 个样本的因变量 y 值，x
(i)
为第 i 个样本的自变量 x 值，由于有 p 个自变量
x
i
，所以 x
(i)
= ( x
(i)
1
, x
(i)
2
, . . . , x
(i)
p
) ∈ R
p
，h
θ
x
(i)
为
h
θ
x
(i)
=
p
j=1
θ
j
x
(i)
j
http://www.ma-xy.com 10 http://www.ma-xy.com
http://www.ma-xy.com
第一章 回归模型 1.2 参数回归
梯度下降算法 GD
Gradient Descent: 在更新每一个参数 θ
j
, j = 1, . . . , p 时，都使用所有样本，则有 θ
j
的梯度
方向
∇θ
j
=
∂J
∂θ
j
= −
1
n
n
i=1
y
(i)
− h
θ
x
(i)
x
(i)
j
于是，参数更新公式为
θ
j
←: θ
j
− η∇θ
j
其中：η 为学习率，可以是变化的学习率 η
t
。由此算法得到的 θ 是一个全局最优解，但是每一步
都要用到所有的样本数据 (训练集)，当 n 很大时，效率非常低。记 ∇
θ
J(θ) 为 J 关于 θ 的梯度。
随机梯度下降算法 SGD
Stochastic Gradient Descent: 每次迭代使用一个样本，并用次样本来更新所有参数 θ
j
(j =
1, . . . , p)，有 θ
j
的梯度方向
∇θ
j
=
y
(i)
− h
θ
x
(i)
x
(i)
j
SGD 伪代码如下 (1)
算法 1 随机梯度下降 SGD
输入: 样本数据 (x
(i)
, y
(i)
)
n
i=1
，x
(i)
∈ R
p
, y
(i)
∈ R，学习率 η
输出: 参数 θ
1: for i = 1, 2, . . . , n do
2: for j = 1, 2, . . . , p do
3: ∇θ
j
=
y
(i)
− h
θ
x
(i)
x
(i)
j
4: θ
j
←: θ
j
− η∇θ
j
5: end for
6: // 或者写为向量形式
7: θ ←: θ −η∇
θ
J(θ; x
i
, y
i
)
8: end for
值得注意的是 SGD 得到的 θ 并不是全局最优解，并且每次 θ
j
下降的方向也不是全局梯度
下降方向，它的下降方向只是使某一样本 x
i
, y
i
的离差平方 (ˆy
i
− y
i
)
2
最小。由于这种异步操作，
SGD 并不易于并行实现，但其收敛速度快，并且，从另一个角度来看，非全局梯度下降有利于跳
出局部极小解，因此，SGD 在非凸优化中也表现突出 (如：神经网络和深度网络等)。
小批量梯度下降算法 MBGD
上面的梯度下降算法 GD 每次迭代使用所有样本，随机梯度下降算法 SGD 每次迭代使用一
个样本，自然地，我们想到一次使用部分样本，我们称这种一次迭代使用部分样本的算法为小批
http://www.ma-xy.com 11 http://www.ma-xy.com
http://www.ma-xy.com
1.2 参数回归 第一章 回归模型
量随机梯度下降算法 (Mini-Batch Gradient Descent)。每次随机的从样本 n 中挑取 b 个样本，这
里的 b 即为批量的大小，用着 b 个样本来更新 θ，有 θ
j
的梯度方向为
∇θ
j
=
1
b
b
i=1
y
(i)
− h
θ
x
(i)
x
(i)
j
MBGD 伪代码如下 (2)
算法 2 小批量梯度下降算法 MBGD
输入: 样本数据 (x
(i)
, y
(i)
)
n
i=1
，x
(i)
∈ R
p
, y
(i)
∈ R，学习率 η，批量大小 b，迭代次数 t，最大迭
代次数 T
输出: 参数 θ
1: for t = 1, 2, . . . , T do
2: 从样本 n 中挑取 b 个样本 x
t
, y
t
3: for j = 1, 2, . . . , p do
4: ∇
θ
j
=
1
b
b
i=1
y
(i)
− h
θ
x
(i)
x
(i)
j
5: θ
j
←: θ
j
− η∇θ
j
6: end for
7: // 或者写为向量形式
8: θ ←: θ −η∇
θ
J(θ; x
t
, y
t
)
9: end for
可以发现，MBGD 是 BGD 和 GD 的综合，相对于 SGD，MBGD 降低了随机性带来的波
动，使更新更加稳健；相对于 BGD，MBGD 提高了收敛速率，同时随机性有助于跳出局部极小
值，因此更多时候，我们采用 MBGD 算法。但如 SGD 那样，MBGD 的更新方向不是全局梯度
方向，其解亦非全局最优解。
从收敛速度来看，SGD 能够在目标函数是强凸且递减学习率的情况下，以 O(
1
T
) 次线性收
敛，而 GD 以 O(ρ
T
)(ρ < 1) 线性收敛，其中：T 为迭代总次数。
动量因子法 Momentum
动量因子法是模拟物理中动量的概念，用累计的动量来替代梯度，可以参考
1999.Qian
3。动
量项为
v
t
= rv
t−1
+ η∇θ
其中：∇θ = (∇θ
1
, . . . , ∇θ
j
, . . . , ∇θ
p
)，动量项超参数 r < 1，一般情况取 0.9 以下。参数更新为
θ ←: θ −v
t
在更新参数时，如果当前 θ
j
的梯度方向与上次梯度方向相同，则 θ
j
在这个方向上的更新会
更快，当达到局部极小值附近时，∇ → 0，r 使得更新幅度增大，从而可以跳出局部解。
http://www.ma-xy.com 12 http://www.ma-xy.com
http://www.ma-xy.com
第一章 回归模型 1.2 参数回归
NAG
1983.Nesterov2提出 Nesterov accelerated gradient，NAG 不仅增加了动量项，并且在计算参
数的梯度时，从损失函数 J(θ) 中减去了动量项
v
t
= rv
t−1
+ η∇
θ
J(θ − rv
t−1
)
θ 的更新公式为
θ ←: θ −v
t
其实，momentum 和 nesterov 项都是为了使梯度更加灵活，但是人工设置学习率还是有些
主观，接下来介绍几种自适应学习算法。
Adagrad
2011.Duchi,Hazan 和 Singer1提出 Adagrad 算法。前面介绍的所有算法，它们在每次更新参
数时，对所有参数 θ
i
, . . . , θ
p
都使用相同的学习率 η，如果数据是稀疏的，或者每个特征 x
j
有着
不同的统计特征，那么，便不能使用相同的学习率。那些很少出现的特征应该使用一个较大的学
习率。
Adagrad 能够对每个参数 θ
j
自适应一个学习率 η
j
t
，对稀疏特征，得到一个大的学习率，因
此，该 Adagrad 算法适合处理有稀疏特征的数据。令参数梯度为
∇θ
j
= ∇
θ
J(θ
j
)
则参数 θ
j
的 Adagrad 更新公式为
θ
j
←: θ
j
− η
1
G
t,j
+ ϵ
∇θ
j
其中：G
t
∈ R
p×p
是一个对角矩阵，其第 j 个对角元素 e
jj
为过去到当前第 j 个参数 θ
j
的梯度
平方和，即
e
jj
←: e
jj
+ ∇θ
2
j
ϵ 是为了使分母不为 0，通常取 1e − 8。另外，如果分母不开根号，算法的性能会非常糟糕。上
式的向量形式为
θ ←: θ −η
1
√
G
t
+ ϵ
· ∇θ
Adadelta
虽然 Adagrad 能够自适应学习率，但其仍然依赖于人工设置的全局学习率 η，并且 Adagrad
还需要计算参数梯度序列平方和 e
jj
，这会使存储量增加，此外，学习率是不断衰减的，最终达
到一个非常小的值。为了克服 Adagrad 的这个缺点，2012.Zeiler6提出了 Adadelta 算法，令
g
t
= ∇
θ
J(θ)
http://www.ma-xy.com 13 http://www.ma-xy.com
http://www.ma-xy.com
1.2 参数回归 第一章 回归模型
其中：t 表示迭代次数。计算梯度平方的期望然后累加
E|g
2
|
t
= E|g
2
|
t−1
+ (1 − γ)g
2
t
计算参数更新量
∆θ
t
=
t−1
r=1
∇θ
r
E|g
2
|
t
+ ϵ
最终，有参数更新公式
θ
t+1
←: θ
t
− ∆θ
t
注意，此时 θ 的更新已经不再依赖于全局学习率了。
RMSprop
RMSprop 是 Adadelte 的中间形式
E|g
2
|
t
= γE|g
2
|
t−1
+ (1 − γ)g
2
t
∆θ
t
= η
1
E|g
2
|
t
+ ϵ
· g
t
θ 的更新公式为
θ
t+1
←: θ
t
− η∆θ
t
在深度学习模型中，Hinton 建议将参数设置为 r = 0.9, η = 0.01。
Adam
Adam(Adaptive Moment Estimation) 利用梯度的一阶矩估计和二阶矩估计动态调整每个参
数的学习率
v
t
= β
1
(v
t−1
) + (1 −β
1
)g
t
n
t
= β
2
(n
t−1
) + (1 −β
2
)g
2
t
ˆv
t
=
v
t
1 − β
t
1
ˆn
t
=
n
t
1 − β
t
2
∆θ
t
= η
ˆv
t
√
ˆn
t
+ ϵ
最终，θ 的更新公式为
θ
t+1
←: θ
t
− ∆θ
t
http://www.ma-xy.com 14 http://www.ma-xy.com
http://www.ma-xy.com
第一章 回归模型 1.3 贝叶斯回归模型
上面介绍了随机梯度下降算法及其变体，我们引入两幅经典的图
(1.4) 来比较随机梯度下降
算法及其变体
(a) (b)
图 1.4: 随机梯度下降算法及其变体的比较
从图 (1.4) 中我们可以看出，在优化问题的鞍点处 (saddle points：某些维度上梯度为 0，某
些不为 0)，SGD、Momentum 和 NAG 一直在鞍点梯度为 0 的方向上震荡，很难打破鞍点位置
的对称性；Adagrad、Adadelta 和 RMSprop 等能够很快的向梯度不为 0 的方向上移动，并且，
这些自适应学习方法有更好的收敛性和收敛速度。
前面，我们提到过 SGD 是次线性收敛的，而 GD 是线性收敛的，那么能否改进 SGD，使
SGD 能够达到线性收敛速度？有的，下面我们来介绍 SAG 和 SVRG。
SAG
SAG 的伪代码如下 (3)
SAG 为每一个样本 (x
i
, y
i
)(伪代码中写为 x
(i)
, y
(i)
) 都维护了一个梯度 g
i
，随机选择一个样
本 x
i
, y
i
来更新 d，并用 d 来更新参数 θ。
SVRG
SVGR 的伪代码如下 (4)
SVRG 的思路是：每隔一段时间计算一次所有样本的梯度 ˜g，每个阶段内部的单次更新采用
∇
θ
J
i
(θ
j−1
) − ∇
θ
J
i
(
˜
θ) + ˜g 来更新参数 θ。SVRG 每次更新最多计算两次梯度，相对 SAG 而言，
不需要在内存中为每个样本维护一个梯度，从而节省了内存。
1.3 贝叶斯回归模型
前面介绍了参数回归的线性回归、参数估计方法、正则回归以及随机梯度下降优化算法，下
面，我们在贝叶斯框架下讨论回归模型的参数估计问题，称基于贝叶斯方法的 (广义) 线性回归
模型为贝叶斯线性回归。
http://sebastianruder.com/optimizing-gradient-descent/?url_type=39&object_type=webpage&pos=1
http://www.ma-xy.com 15 http://www.ma-xy.com
http://www.ma-xy.com
1.3 贝叶斯回归模型 第一章 回归模型
算法 3 Basic SAG method for minimizing J(θ) =
1
n
n
i=1
f
i
(θ)with step size η
输入: 样本数据 (x
(i)
, y
(i)
)
n
i=1
，x
(i)
∈ R
p
, y
(i)
∈ R，学习率 η，迭代次数 t，最大迭代次数 T
输出: 参数 θ
1: 初始化。d = 0, g
i
= 0
2: for t = 1, 2, . . . , T do
3: 从样本 n 中挑取 1 个样本 x
(i)
, y
(i)
4: for j = 1, 2, . . . , p do
5: d
j
= d
j
− g
j
i
+ ∇
θ
J(θ
j
|x
(i)
, y
(i)
)
6: g
j
i
= ∇
θ
J(θ
j
|x
(i)
, y
(i)
)
7: θ
j
←: θ
j
−
η
n
d
j
8: end for
9: // 或者写为向量形式
10: d = d − g
i
+ ∇
θ
J
i
(θ)
11: g
i
= ∇
θ
J
i
(θ)
12: θ ←: θ −
η
n
d
13: end for
1.3.1 贝叶斯统计基础
在概率论中，一直存在两大学派：频率学派和贝叶斯学派。就参数估计问题而言，频率学派
认为参数 θ 是未知待求的固定量，而贝叶斯学派认为参数 θ 是随机变量，因此，Bayes 的参数会
存在分布。称在没给出样本前得到的分布为先验分布，先验分布是由经验积累得到的信息，如：
一个企业家认为“一项新产品在未来市场上的畅销”的概率为 0.8，这里的 0.8 就是企业家根据
自己多年的经验得到的结论，即先验概率，但是后期是否真的为 0.8，还需要结合实际抽样来下
结论。
为了便于叙述贝叶斯参数估计，我们先来简单介绍一下贝叶斯公式。设 A, B 和 A
i
是事件，
P (A) 为事件 A 发生的概率，P (B|A
i
) 表示 A
i
发生导致 B 发生的概率，并假设可以导致 B 发
生的事件有 A
1
, A
2
, . . . , A
i
, . . . ，有如下条件概率公式
P (A|B) =
P (AB)
P (B)
P (A
i
B) = P (B)P (A
i
|B) = P (A
i
)P (B|A
i
)
由此，可以得到离散形式的贝叶斯公式
P (A
i
|B) =
P (A
i
)P (B|A
i
)
P (B)
=
P (A
i
)P (B|A
i
)
i
P (A
i
)P (B|A
i
)
(1.9)
通俗的说，贝叶斯公式的概率意义是：事件 B 发生了，事件 B 发生是由事件 A
i
引起的概
率等于 A
i
发生且 A
i
发生导致 B 发生的概率除上所有可能导致 B 发生的情况。
http://www.ma-xy.com 16 http://www.ma-xy.com
http://www.ma-xy.com
第一章 回归模型 1.3 贝叶斯回归模型
算法 4 Basic SVRG method for minimizing J(θ) =
1
n
n
i=1
f
i
(θ)with step size η
输入: 样本数据 (x
(i)
, y
(i)
)
n
i=1
，x
(i)
∈ R
p
, y
(i)
∈ R，学习率 η，迭代次数 t，最大迭代次数 T
输出: 参数 θ
初始化：
˜
θ
0
for t = 1, 2, . . . , T do
˜
θ ←
˜
θ
t−1
˜g =
1
n
n
i=1
∇f
i
(
˜
θ)
从样本 n 中挑取 1 个样本 x
(i)
, y
(i)
for j = 1, 2, . . . , p，p ≡ m do
d
j
= d
j
− g
j
i
+ ∇
θ
J(θ
j
|x
(i)
, y
(i)
)
∆θ
j
= η
∇
θ
J
i
(θ
j−1
) − ∇
θ
J
i
(
˜
θ) + ˜g
θ
j
←: θ
j
− ∆θ
j
end for
option 1: set
˜
θ
t
=
˜
θ
p
option 2: set
˜
θ
t
=
˜
θ
j
，for randomly chose j ∈ {0, 1, . . . , p −1}
end for
贝叶斯学派认为未知待求参数是随机的，有一定的分布。我们令 θ 为未知参数，记其先验概
率为 π(θ)，令 x ∈ R
p
为变量，x = (x
1
, . . . , x
n
) 为样本，则样本 x 的联合概率密度函数是在 θ
给定下的条件密度，故 x 的密度函数为 p(x|θ)，这里将 θ 视为导致样本 x 出现的原因，相当于
A。由上面的贝叶斯公式 (1.9)，我们有 θ 的后验分布
π(θ|x) =
p(x|θ)π(θ)
Θ
p(x|θ)π(θ)dθ
=
p(x|θ)π(θ)
p(x)
(1.10)
其中：π(θ|x) 就是在样本 x 给定后，θ 的后验概率，我们要求使得 π(θ|x) 最大的 θ
∗
，这样的参
数估计方法称为最大后验概率方法。
max
θ
π(θ|x)
由于式 (1.10) 分母中的样本联合概率分布 p(x) 相对于 θ 独立，故
max
θ
π(θ|x) = max
θ
p(x|θ)π(θ)
注意：这里的 p(x|θ) 和极大似然估计中的 p(x; θ) 是一样。一般而言，我们并没有任何辅助
信息来确定先验概率密度 π(θ)，所以 π(θ) 是主观意义上的，可以取正态分布，Gamma 分布族
等等。思考：如果按照贝叶斯统计的思想，我们给定参数 θ 的带参数先验分布 π(θ|α, β)，那么
α, β 也是随机变量，也会有相应的先验分布和后验分布，如此下去，以致无穷，我们称 α, β 为超
参数，后面会介绍这种“无穷”问题的处理方法。
http://www.ma-xy.com 17 http://www.ma-xy.com
http://www.ma-xy.com
1.3 贝叶斯回归模型 第一章 回归模型
1.3.2 贝叶斯线性回归模型
在前面的线性回归模型中，我们假设 ε
i
iid
∼ N(0, σ
2
)，y
i
|x
i
iid
∼ N(ˆy
i
, σ
2
)。并且，已经有了在
参数 θ = (w, σ
2
), w ∈ R
m
, σ ∈ R 给定下样本的联合概率密度
p(y|x, w, β) =
n
i=1
N(w
T
ϕ(x
i
), β
−1
)
其中：β = (σ
2
)
−1
。此联合概率分布可以记为 p(y|θ)，即前面说的 p(x|θ)，那么，根据贝叶斯公
式 (1.9)，θ 的后验概率 π(θ|y) 为
π(θ|y) =
p(y|θ)π(θ)
Θ
p(y|θ)π(θ)dθ
我们的目标是求 θ，使得后验概率 π(θ|y) 最大，但还没有给出先验分布 π(θ)。下面就来讨论
π(θ)。给出参数 w 的先验概率分布，在这期间，把 β 当作已知常数，假设 1：w = (w
1
, w
2
, . . . , w
p
)
的 w
i
是独立同分布的 (p ≡ m)，或者说它们的先验分布是一样的；假设 2：w
i
∼ N(0, α
−1
)，其
中：α 为超参数。
由各参数 w
i
先验分布一样，为正态分布 N(0, α
−1
)，我们有
p(w
i
; α) = N (0, α
−1
)
π(w) = π(w; α) = N (0, α
−1
I)
这里的 π(w; α) 也可以记为 π(w|α)，表示分布 π(w) 是由参数 α 控制的。给出了先验 π(w|α) 之
后，还要计算样本条件密度函数 p(y|θ) 以及后验概率密度 π(θ|y)。p(y|θ) 即为样本的条件联合概
率密度，而对于 π(θ|y) 的求解，更详细的推导可以参考：《高等数理统计 (第二版)》茆诗松 P350
例
5.26
或者《
PRML
》
P67-P69
。下面，我们用一示例来展示后验分布
π
(
θ
|
y
)
为正态分布。
示例 (正态均值的共轭先验分布为正态分布) 我们用此示例来展示：在方差已知的情况下，正
态均值的共轭先验分布为正态分布。设 y = (y
1
, y
2
, . . . , y
n
) 是来自正态分布 N (θ , β
−1
) 的样本
(已知，y 是来自 N (w
T
ϕ(x), β
−1
) 的样本，而 w 就是 θ)，其中 β
−1
= σ
2
已知，且设参数 (均值)θ
的先验分布为正态分布 θ ∼ N(µ, τ
2
)，即
π(θ) =
1
√
2πτ
e
−
(θ−µ)
2
2τ
2
其中：µ, τ 已知。由此，可以写出样本 y 与 θ 的联合密度函数
h(y, θ) = k
1
exp
−
1
2
nθ
2
− 2nθ¯y +
n
i=1
y
2
i
σ
2
+
θ
2
− 2µθ + µ
2
τ
2
其中：k
1
= (2 π)
−
n+1
2
τ
−1
σ
−n
，¯y = E(y)。再记
σ
2
0
=
σ
2
n
, A =
1
σ
2
0
+
1
τ
2
, B =
¯y
σ
2
0
+
µ
τ
2
, C =
n
i=1
y
2
i
σ
2
+
µ
2
τ
2
http://www.ma-xy.com 18 http://www.ma-xy.com
http://www.ma-xy.com
第一章 回归模型 1.3 贝叶斯回归模型
则有
h(y, θ) = k
1
−
1
2
[Aθ
2
− 2θB + c]
= k
1
exp
−
(θ − B/A)
2
2/A
−
1
2
(C − B
2
/A)
(1.11)
容易计算样本 y 的边缘分布为
m(y) =
∞
−∞
h(y, θ)dθ
= k
1
exp
−
C − B
2
/A
2
2π
A
1
2
(1.12)
上式 (1.11) 和 (1.12) 相除，可以得到 θ 的后验分布
π(θ|y) =
2π
A
−
1
2
exp
−
(θ − B/A)
2
2/A
即后验分布 π(θ|y) 为正态分布，记
π(θ|y) = N (µ
1
, σ
2
1
)
其中：µ
1
=
B
A
，σ
2
1
= ( σ
2
0
+ τ
−2
)
−1
。
我们不加证明的给出上例中的向量形式的结论：先给出 θ 的先验分布和 y 的联合概率密度
π(θ) = N (µ, Λ
−1
)
p(y|θ) = N (ˆy = Ax + b, L
−1
)
则 y 的边缘分布 m(y) 以及 θ 的后验分布 π(θ|y) 为
m(y) = N (Aµ + b, L
−1
+ AΛ
−1
A
T
)
π(θ|y) = N (Σ(A
T
L(y − b) + Aµ), Σ)
其中：Σ = (Λ + A
T
LA)
−1
。
上面，为了便于区分 y 和 θ，我们将 θ 的先验分布和后验分布写为 π(θ) 和 π(θ|y)，将 y 的
条件联合概率密度写为 p(y|θ) = p(y; θ)，y 的概率分布 (边缘分布) 写为 m(y)。我们将 π(θ|y) 记
回 p(θ|y)，经过上面的介绍，对于参数 w 的后验分布 π(w|y)，我们有
p(w|y) = N (m
N
, S
N
)
其中：N ≡ n 为样本数，
m
N
= βS
N
ϕ
T
y
S
−1
N
= αI + βϕ
T
ϕ
http://www.ma-xy.com 19 http://www.ma-xy.com
http://www.ma-xy.com
1.3 贝叶斯回归模型 第一章 回归模型
我们的目标是求使得后验概率最大的 θ，那么在 β
−1
给定的情况下，可以求 w，使得 p(w|y)
最大，由于 p(w|y) 是正态的，所以和极大似然估计相似
max
w
L(w) = p(w|y)
对上式取对数，极大情况不变，有
max
w
L(w) ≡ max
w
ln L(w)
= max
w
ln p(w|y)
= max
w
−
β
2
n
i=1
y
i
− w
T
ϕ(x
i
)
2
−
α
2
w
T
w + constant
其中：constant 是与 w 无关的常数。如果我们已知 α, β 超参数，并且根据上面的优化模型求出
了最优的
w
。那么，接下来的问题是：给定一个新的样本
x
∗
，如何预测 y
∗
，即
p(y
∗
|y, x, α, β) =
w
p(y
∗
|x
∗
, w, β)p(w|y, α, β)dw
注意到上式的预测是两个高斯分布的卷积，p(w|y, α, β) 为后验概率，前面已经介绍了，p(y
∗
|x
∗
, w, β)
为条件概率分布。由卷积公式，可以求得预测分布的形式
p(y
∗
|x, y, α, β) = N(y|m
T
N
, σ
2
N
(x))
其中：σ
2
N
= β
−1
+ ϕ(x)
T
S
N
ϕ(x)，σ
2
N
为预测分布的方差，其第二项 ϕ(x)
T
S
N
ϕ(x) 反应了参数 w
关联的不确定性。
由于噪声 ε 和 w 的分布是相互独立的高斯分布，因此，其值有
σ
2
N+1
(x) ⩽ σ
2
N
(x)
且当样本数 N → ∞ 时，σ
2
N
(x) 只与 β 有关，ϕ(x)
T
S
N
ϕ(x)。
上面在计算预测分布时，假设已知超参数 α, β，但在实际中，大多数情况下并不知道引入的
超参数 α, β 是多少，这使得问题难以解决。当然，可以像 w, σ 那样，对超参数再引入超先验分
布，那样，预测分布可以通过对 w, α, β 求积得到
p(y
∗
|y) =
p(y
∗
|w, β)p(w|y, α, β)p(α, β|y)dwdαdβ
并且，如果定义 α 和 β 上的共轭先验分布为 Gamma，那么，上式可以解析计算出来，得到
w 上的学生 t 分布。下面介绍一种求解 α, β 的方法 - 最大边缘似然函数法。
最大边缘似然函数法
最大边缘似然函数法又称为第二类最大似然方法 (广义极大似然或者证据近似)。首先，给出
边缘必然函数 p(y|α, β) 的计算
p(y|α, β) =
w
p(y|w, β)p(w|α)dw
上式是一个高斯分布与高斯分布的卷积，可以按下面的线性高斯模型的条件概率分布求解。
http://www.ma-xy.com 20 http://www.ma-xy.com
http://www.ma-xy.com
第一章 回归模型 1.3 贝叶斯回归模型
引理 (高斯分布的卷积公式) 由
p(w) = N (0, Λ
−1
)
p(y|w) = N (Ax + b, L
−1
)
有
p(y) = p(w)p(y|w) = N (Aµ + b, L
−1
+ AΛ
−1
A
T
)
根据上面的高斯分布卷积公式，我们将 p(y|α, β) 显式的写出来
p(y|α, β) =
β
2π
n
2
α
2π
n
2
exp{E(w)}dw
其中：
E(w) =
β
2
||y − ϕw||
2
+
α
2
w
T
w
=
β
2π
n
2
α
2π
n
2
exp{−E(m
N
)}(2π)
M
2
|A|
−
1
2
其中：
E(m
N
) =
β
2
||y − ϕm
N
||
2
+
β
2
m
T
N
m
N
m
N
= βA
−1
ϕ
T
y
A = S
−1
N
= αI + βϕ
T
ϕ
这里的 N ≡ n 为样本数，M ≡ m 为映射 ϕ 的个数，ϕ = (ϕ
1
, ϕ
2
, . . . , ϕ
m
)。
我们的目标是：求 α, β，使边缘似然函数 p(y|α, β) 最大。由于 p(y|α, β) 为高斯分布，那么
和极大似然估计相似，对其取对数，有
max
α,β
ln L(α, β) = ln p(y|α, β)
=
M
2
ln α +
N
2
ln β −E(m
N
) −
1
2
ln |A| −
N
2
ln(2π) (1.13)
对于上面的边缘似然目标 (1.13)，我们先考虑 p(y|α, β) 关于 α 的最大化，然后再考虑其关
于 β 的最大化。
(1) 先来考虑 p(y|α, β) 关于 α 的最大化。由 A = αI + βϕ
T
ϕ 可知，A 的特征值为 α + λ
i
。
现在考虑式 (1.13) 中 ln |A| 关于 α 的导数
d ln |A|
dα
=
d
dα
ln
i
(λ
i
+ α)
=
d
d
α
i
ln(λ
i
+ α)
=
i
1
λ
i
+ α
http://www.ma-xy.com 21 http://www.ma-xy.com
http://www.ma-xy.com
1.3 贝叶斯回归模型 第一章 回归模型
所有
∂ ln L(α, β)
∂α
=
M
2α
−
1
2
m
T
N
m
N
−
1
2
i
1
λ
i
+ α
令上式的导数为 0，可以求解 α
M
2α
−
1
2
m
T
N
m
N
−
1
2
i
1
λ
i
+ α
= 0
上式两边乘 2α，有
αm
T
N
m
N
= M − α
i
1
λ
i
+ α
= r
由于 i 的求和项共有 M 项
M =
i
λ
i
+ α
λ
i
+ α
所以，有
r =
i
λ
i
λ
i
+ α
α =
r
m
T
N
m
N
注意到这是一个求解 α 的隐式解，不仅 r 与 α 相关，后验概率本身的众数 m
N
也与 α 有
关，因此使用迭代方法来求解 α。
Step1. 设定 α 的初始值 α
0
。
Step2. 求 m
N
。
m
N
= βS
N
ϕ
T
y
S
−1
N
= αI + βϕ
T
ϕ
A = S
−1
N
Step3. 计算 r。
r =
i
λ
i
λ
i
+ α
其中：α + λ
i
为 A 的特征值。
Step4. 求解 α。
α ←:
r
m
T
N
m
N
注意：由于矩阵 ϕ
T
ϕ 是固定的，因此，在最开始计算一次特征值后，接下来只需要乘以 β 就可
以得到 λ
i
的值。
http://www.ma-xy.com 22 http://www.ma-xy.com
http://www.ma-xy.com
第一章 回归模型 1.4 RVM 稀疏向量机
(2) 上面，给出了 p(y|α, β) 关于 α 的最大化。下面，给出 p(y|α, β) 关于 β 的最大化。注意
到特征值 λ
i
正比于 β，因此
d
dβ
=
λ
i
β
于是
d
dβ
ln |A| =
d
β
i
ln(λ
i
+ α)
=
1
β
i
λ
i
λ
i
+ α
=
r
β
继而，式 (1.13) 关于 β 的导数为
∂ ln L(α, β)
∂β
=
N
2β
−
1
2
n
i=1
y
i
− m
T
N
ϕ(x
i
)
2
−
r
2β
令上面的导数为 0，整理后，有
1
β
=
1
N −r
n
i=1
y
i
− m
T
N
ϕ(x
i
)
2
与 α 一样，β 的计算公式也是一个隐式的，其迭代算法表述为：
Step1.β
0
。
Step2. 计算 m
N
, r。
Step3. 计算 β。
1.4 RVM 稀疏向量机
1.4.1 RVM 的建立
2001.Tipping4提出了一个用于回归和分类问题的稀疏贝叶斯方法 - RVM?。前面我们构建的
广义线性回归模型为
y =
M−1
i=0
w
i
ϕ
i
(x) + ε
向量形式为
y = w
T
ϕ(x) + ε
其中：ε 为误差项，ϕ(x) 为设计矩阵，ϕ(x) = (ϕ
0
(x), ϕ
1
(x), . . . , ϕ
M
−
1
(x))
T
，ϕ
i
(x) 是非线性基函
数，M ≡ m 为基函数个数，N ≡ n 为样本个数。
http://www.ma-xy.com 23 http://www.ma-xy.com
http://www.ma-xy.com
1.4 RVM 稀疏向量机 第一章 回归模型
特别地，当基函数 ϕ(x) 由核函数给出，样本中的每一个样本点都会有一个核，即基函数个
数 M = n，则上面的广义线性回归模型可以写为如下形式
y =
n
i=1
w
i
K(x, x
i
) + w
0
+ ε
上述模型参数数量为 n + 1，不包含 ε 的方差 σ
2
。同前面的假设一样，我们假设 ε
iid
∼ N(0, σ
2
)，
令 β = (σ
2
)
−1
，则由此可以写出样本 x, y 的似然函数
p(y|x, w, β) =
n
i=0
p(y
i
|x
i
, w, β)
在贝叶斯线性回归当中，我们假设 ∀i, j ∈ [0, n]w
i
, w
j
有相同的分布 w
i
∼ N(0, α
−1
)，或者整体
写为 p(w|α) = N(0, α
−1
I)。现在，我们不要求 w
i
, w
j
具有相同的高斯先验分布，假设 w
i
仍然是
高斯分布，但是它们的方差不同
p(w
i
|α
i
) = N(0, α
−1
i
) i = 0, 1, . . . , n
并且假设随机变量 w
i
, w
j
之间相互独立，同样我们可以写出 w 的分布
p(w|α) =
n
i=0
N(0, α
−1
i
)
其中：α
i
表示随机变量 w
i
的精度，即方差倒数；α = (α
0
, α
1
, . . . , α
n
)
T
上面给出了参数 w 的先验分布，由前面的知识，我们知道参数 w 的后验分布是一个高斯分
布
p(w|y, x, α, β) = N(µ, Σ)
=
p(y|wα, β)p(w, α, β)
p(y)
=
p(y|wα, β)p(w, α, β)
p(y|w, α, β)p(w, α, β)dwdαdβ
其中：
µ = βΣϕ
T
y
Σ = (A + βϕ
T
ϕ)
−1
A = diag(a
i
)
ϕ 为 N ×N 设计矩阵。目标是求使得后验概率 p(w|y, x, α, β) 最大的 w ，但是并不知道方差倒数
β 和超参数 α 的信息，下面来求 α, β。由于是在 Bayes 框架下，参数 w 是随机的，我们也假设
α, β 是随机的，且有一定的分布。由于 Gauss 分布方差倒数 β 的共轭先验分布为 Gamma 分布，
因此，假设 α
i
, σ
2
的超先验分布为
p(α) =
N
i=0
Gamma(α|a, b)
p(β) = Gamma(c, d)
http://www.ma-xy.com 24 http://www.ma-xy.com
http://www.ma-xy.com
第一章 回归模型 1.4 RVM 稀疏向量机
其中：
Gamma(a, b) = Γ(a)
−1
b
a
α
a−1
e
−bα
Γ(a) =
∞
0
t
a−1
e
−t
dt
1.4.2 最大边缘似然方法
仿照前面求解 w 的最大后验概率方法 (最大边缘似然方法) 来求解超参数 α, β。超参数 α, β
的后验分布为
p(α, β|y) ∝ p(y|α, β)p(α)p(β) (1.14)
则可以求 α, β，使 α, β 的后验概率 (1.14) 最大。对于式 (1.14)，先来计算 p(y|α, β)，p(y|α, β)
即为前面所说的边缘似然函数
p(y|α, β) =
p(y|x, w, β)p(w|α)dw
= N(µ
T
ϕ(x), Σ)
= (2 π)
−
n
2
|β
−1
I + ϕA
−1
ϕ
T
|
1
2
exp
−
1
2
y
T
(β
−1
I + ϕA
−1
ϕ
T
)
−1
y
其中：Σ = β
−1
I + ϕ
T
A
−1
ϕ，A = diag(a
0
, a
1
, . . . , a
n
)。所以求 α, β 使 p(y |α, β) 最大，可以写为
max
α,β
p(y|α, β) ≡ max
α,β
p(y|α, β)p(α)p(β) (1.15)
对上面的目标函数 (1.15) 取对数 log，对于 α, β，忽略本身对数独立性，并且为了便利，我们在
log α, log β 上求最大化。因为 p(log α) = αp(α)，于是目标函数 (1.15) 变为
L(α, β) = log p(y|log α, log β) +
N
i=0
log p(log α
i
) + log p(log β)
= −
1
2
log |β
−1
I + ϕA
−1
ϕ
T
|y
T
(β
−1
I + ϕA
−1
ϕ
T
)
−1
+
N
i=0
(α log α
i
− bα
i
) + c log β − dβ (1.16)
在编程中，一般将 a, b, c, d 设置为 10
−4
。我们求 α, β 使式 (1.16)
L
(
α, β
)
最大，并且在
log
α,
log
β
上求最大化，求导并令导数为 0，有
∂L
∂ log α
i
= 0
∂L
∂ log β
= 0
为了求解上面的等式，我们对式 (1.16) 进行逐项分析。令 C = β
−1
I + ϕA
−1
ϕ
T
(1)
log |C| = log |β
−1
I + ϕA
−1
ϕ
T
|
= −log |Σ| − N log β − log |A|
http://www.ma-xy.com 25 http://www.ma-xy.com
http://www.ma-xy.com
1.4 RVM 稀疏向量机 第一章 回归模型
其中：Σ = (A + βϕ
T
ϕ)
−1
是 N ×N 矩阵。
(2) 由 Woodbury 倒置恒等式：
|A||β
−1
I + ϕA
−1
ϕ
T
| = |β
−1
I||A + β
T
ϕ|
有
(β
−1
I + ϕA
−1
ϕ
T
)
−1
= βI − βϕ(A + βϕ
T
ϕ)
−1
ϕ
T
β
于是
y
T
(β
−1
I + ϕA
−1
ϕ
T
)y = βy
T
y − βy
T
ϕΣϕ
T
βy
= βy
T
(y − ϕµ)
其中：µ 在前面已经说过了，µ = βΣϕ
T
y 是后验概率 p(w|y, α, β) 的均值向量。
βy
T
(y − ϕµ) = β||y − ϕµ||
2
+ βy
T
ϕµ − βµ
T
ϕ
T
ϕµ
= β||y − ϕµ||
2
+ µ
T
Σ
−
1
µ − βµ
T
ϕ
T
ϕµ
= β||y − ϕµ||
2
+ µ
−1
Aµ
注：Σ = (A + βϕ
T
ϕ)
−1
的计算复杂度为 O(N
3
)。经过上面的分析，现在，我们可以给出求导的
结果：
¬
∂L
∂ log α
i
=
1
2
[1 − α
i
(µ
2
i
+ Σ
ii
)] + a + bα
i
令上式为 0，可以得到 α
i
的更新迭代公式
α
i
=
1 + 2a
µ
2
i
+ Σ
ii
+ 2b
令 r
i
= 1 − α
i
Σ
ii
，则
α
i
=
r
i
+ 2a
µ
2
i
+ 2b
其中：µ
i
视为后验概率分布 p(w|y, α, β) 的均值向量 µ 的第 i 个分量；Σ
ii
为后验概率 p(w|y, α, β)
的方差 Σ 的第 i 个对角元素。
∂L
∂ log β
=
1
2
N
β
− ||y − ϕµ||
2
− tr(Σϕ
T
ϕ)
+ c − dβ
这里，tr(Σϕ
T
ϕ) 也可以写为 β
−1
Σ
i
r
i
。令上式为 0，可以得到 β 的更新迭代公式为
β
−1
=
||y − ϕµ||
2
+ 2d
N −Σ
i
r
i
+ 2c
http://www.ma-xy.com 26 http://www.ma-xy.com
http://www.ma-xy.com
第一章 回归模型 1.4 RVM 稀疏向量机
注：在迭代过程中，许多 α
i
→ ∞，相应的权重 w
i
→ 0，与之相应的基函数 φ
i
≡ k
i
不起作用，
被删除，从而实现稀疏化。
设通过上面的迭代算法求解得到的最优超参数为 α
∗
, β
∗
，则对于一个新的样本 x
∗
，我们可以
给出 y
∗
的预测分布
p(y
∗
|x
∗
, y, α
∗
, β
∗
) =
p(y
∗
|x
∗
, w, β
∗
)p(w|x, y, α
∗
, β
∗
)dw
= N(µ
T
ϕ(x), σ
2
(x))
其中：
σ
2
(x) = (β
∗
)
−1
+ ϕ(x)
T
Σϕ(x)
Σ = (A + β
∗
ϕTϕ)
−1
1.4.3 EM 算法
上面是基于最大边缘似然来求解 α, β 的，下面，介绍用 EM 算法来求解 α, β。在上面，我
们的目标是求 α, β，使得后验概率 (1.14)p(α, β|y) 或者 p(y|α, β)p(α)p(β) 最大。由于 w 已经被
积分除去了，我们将 w 视为一个潜变量，利用 EM 算法来最优化目标：
¬在 E 步，计算当前 α, β 下 w 的后验概率分布 p(w|y, α, β)，然后找到完整数据下对数似然函数
的期望。
在 M 步，关于 α, β 最大化上述期望。
具体的，对
RVM
模型而言
E[p(w|y, α, β)] = E
w|y,α,β
[log p(y|w, β)p(w|α)p(α)p(β)]
对于 α，忽略 α 本身对数独立性，在 M 步，有
max
α,β
L(α, β) = E[p(w|y, α, β)] = E
w|y,α,β
[log p(w|α)p(α)]
(1) 对于 α。忽略 α 本身对数独立性，等价目标为
max
α
L(α) = E
w|y,α,β
[log p(w|α)p(α)]
将上述 L(α) 关于 α 求导，并令导数为 0，得到 α 的更新公式
α
i
=
1 + 2α
w
2
i
+ 2b
其中：
w
2
i
≡ E
w|y,α,β
[w
2
i
] = Σ
ii
+ µ
2
i
。
(2) 对于 β。忽略其本身的对数独立性，等价目标为
max
β
L(β) = E
w|y,α,β
[log p(y|w, β)p(β)]
将上面的 L(β) 关于 β 求导，并令导数为 0，得到 β 的更新公式为
β
−1
=
||y − ϕµ||
2
+ β
−1
Σ
i
r
i
+ 2d
N + 2c
http://www.ma-xy.com 27 http://www.ma-xy.com
http://www.ma-xy.com
1.4 RVM 稀疏向量机 第一章 回归模型
1.4.4 快速序列算法
2003.Tipping5给出了一种求解 α, β 的快速序列算法。前面我们提到的两种求解参数 α, β 方
法 (最大边缘似然/证据近似和 EM 算法) 在参数求解的过程中，都涉及到一个 N ×N 矩阵求逆
的过程
Σ = (A + βϕ
T
ϕ)
−1
其中：N 为基函数个数 (样本个数，每个样本都有一个基函数)。该矩阵求逆的计算复杂度为
O(N
3
)，存储空间为 O(N
3
)。为了使 RVM 能够快速的学习，Tipping 于 2003 年提出了快速序
列算法，显式的写出边缘似然函数中所有对特定的 α
i
的依赖关系，然后显式的确定驻点。在前
面的 RVM 模型中，假设每一个样本都有一个核函数，从而核函数的个数为 N。现在，我们将核
(基) 函数的个数一般化，设基函数个数为 M，有
y = ϕw + ε
其中：ϕ = (ϕ
1
, ϕ
2
, . . . , ϕ
M
) 是 N ×M 的设计矩阵。因为 ϕ 是 N ×M 的，所以 Σ = (A+βϕ
T
ϕ)
−1
是 M ×M 的矩阵求逆。假设 ε
iid
∼ N (0, σ
2
)，样本的概率分布为
p(y|w, β) = (2π)
−
N
2
σ
−N
exp
−
||y − ˆy||
2
2σ
2
其中：ˆy = ϕw。假设超参数 α
i
, α
j
相互独立，则 p(w) 的先验分布为
p(w|α) = (2π)
−
M
2
M
m=1
α
1
2
m
exp
−
α
m
w
2
m
2
其中；α = (α
1
, α
2
, . . . , α
M
)
T
。由贝叶斯公式，得到 p(w) 的后验分布 p(w|y)
p(w|y, α) = N (µ, Σ)
其中：Σ = (A + βϕ
T
ϕ)
−1
，µ = βΣϕ
T
y，A = diag(α
1
, α
2
, . . . , α
M
)。
稀疏贝叶斯模型中的 α, β 可以用第二类极大似然方法求解，它的目标是求 α, β，使得它们
的边缘似然函数最大。同时，由于本身的最大化与最大化 log 等价，故有
max
α
L(α) = log p(y|α, β)
= log
∞
−∞
p(y|w, β)p(w|α)dw
= −
1
2
[N log 2π + log |C| + y
T
C
−1
y]
其中：C = β
−1
I + ϕA
−1
ϕ
T
。我们将 C 重新写为单一 α
i
, i = 1, 2, . . . , M 的形式，有
C = β
−1
I +
m
α
−1
m
ϕ
m
ϕ
T
m
= β
−1
I +
m=i
α
−1
m
ϕ
m
ϕ
T
m
+ α
−1
i
ϕ
i
ϕ
T
i
= C
−i
+ α
−1
i
ϕ
i
ϕ
T
i
http://www.ma-xy.com 28 http://www.ma-xy.com
http://www.ma-xy.com
第一章 回归模型 1.4 RVM 稀疏向量机
其中：C
−i
是 C 去除第 i 个基函数影响后的矩阵。并且有
|C| = |C
−i
| |1 + α
−1
i
ϕ
T
i
C
−1
−i
ϕ
i
|
C
−1
= C
−1
−i
−
C
−1
−i
ϕ
i
ϕ
T
i
C
−1
i
α
i
+ ϕ
T
i
C
−1
−i
ϕ
i
我们重新写 L(α)，有
L(α) = −
1
2
[N log 2π + log |C| + y
T
C
−1
y]
= −
1
2
N log 2π + log |C
−i
| + y
T
C
−1
−i
y − log α
i
+ log(α
i
+ ϕ
T
i
C
−1
−i
ϕ
i
) −
(ϕ
T
i
C
−1
−i
y)
2
α
i
+ ϕ
T
i
C
−1
−i
ϕ
i
= L(α
−i
) +
1
2
log α
i
− log(α
i
+ ϕ
T
i
C
−1
−i
ϕ
i
) +
(ϕ
T
i
C
−1
−i
y)
2
α
i
+ ϕ
T
i
C
−1
−i
ϕ
i
= L(α
−i
) + ℓ(α
i
)
其中：L(α
−i
) 是从模型中去除 α
i
(和 w
i
, ϕ
i
) 的对数边际似然函数，我们现在将函数 ℓ(α
i
) 中的
项 α
i
分离出来。
求 α 使得 L(α) 最大。L(α) 关于 α 求导，令导数为 0，从而得到 α 的更新公式。
∂L(α)
∂α
i
=
dℓ(α
i
)
dα
i
=
1
2
1
α
i
−
1
α
i
+ ϕ
T
i
C
−1
−i
ϕ
i
−
(ϕ
T
i
C
−1
−i
y)
−2
(α
i
+ ϕ
T
i
C
−1
−i
ϕ
i
)
2
=
α
−1
−i
s
2
i
− (q
2
i
− s
i
)
2(α
i
+ s
i
)
2
其中：q
i
= ϕ
T
i
C
−1
−i
y，s
i
= ϕ
T
i
C
−1
−i
ϕ
i
。令
∂L(α)
∂α
i
= 0 ，求解 α
i
：
1. 如果 q
2
i
> s
i
，则 α
i
=
s
2
i
q
2
i
−s
i
；
2. 如果 q
2
i
< s
i
，则 α
i
→ ∞。(这种情况在前面已经说明了)
注：2002.Tipping 分析了 L(α) 的二阶导数，确定了这些解的唯一性。通过上述方法，可以直接
计算出所有基函数 ϕ
i
的 s
i
, q
2
i
。如果设
S
i
= ϕ
T
i
C
−1
ϕ
i
Q
i
= ϕ
T
i
C
−1
y
则有
s
i
=
α
i
S
i
α
i
− S
i
q
i
=
α
i
Q
i
α
i
− Q
i
http://www.ma-xy.com 29 http://www.ma-xy.com
http://www.ma-xy.com
1.4 RVM 稀疏向量机 第一章 回归模型
当 α
i
→ ∞ 时，s
i
= S
i
, q
i
= Q
i
。并且，在实际的学习过程中，可以利用 Woodbury 恒等式来求
解 S
i
, Q
i
S
i
= ϕ
T
i
Bϕ
i
− ϕ
T
i
BϕΣϕ
T
Bϕ
i
Q
i
= ϕ
T
i
By − ϕ
T
i
BϕΣϕ
T
By
其中：B = β
−1
I。
注：1. 对分类问题，B = diag(β
1
, β
2
, . . . , β
N
)，y = ϕw
MP
+ B
−1
(ˆy −y)；2.Woodbury 恒等式为
(I + AB)
−1
= A(I + BA)
−1
(A + BD
−1
C)
−1
= A
−1
− A
−1
BCD + (A
−1
B)
−1
CA
−1
综上，对一个 RVM 模型，其快速序列算法求解 α, β 的过程为
Step1. 初始化。初始 β 或者 σ
2
。
Step2. 初始化 1 个基函数 ϕ
i
，指定其超参数为 α
i
α
i
=
||ϕ
i
||
2
||ϕ
T
i
y||
2
/||ϕ||
2
− β
−1
其余超参数 α
j
设置为 0。
Step3. 计算 Σ 和 µ，同时计算出所有基函数 ϕ
m
对应的 s
m
, q
m
(for all M base ϕ
m
)。
Step4. 选择候选基函数 ϕ
i
。
Step5. 计算 θ
i
= q
2
i
− s
i
。
Step6. 如果 θ > 0 并且 α
i
< ∞，则模型中基函数 ϕ
i
已存在，并更新 α
i
a
i
=
s
2
i
θ
i
Step7. 如果 θ
i
> 0 并且 α
i
= ∞，则增加 ϕ
i
到模型的基函数中，并更新
a
i
=
s
2
i
θ
i
Step8. 如果 θ ⩽ 0 并且 α
i
< ∞，则在模型中删除基函数 ϕ
i
，并更新 α
i
= ∞。
Step9. 对于回归问题，更新测量误差的方差
σ
2
=
||ˆy − y||
2
N −M + Σ
ii
α
i
Σ
ii
Step10. 更新 Σ, µ 和所有的 s
m
, q
m
s
m
=
α
m
S
m
α
m
− S
m
q
m
=
α
m
Q
m
α
m
− Q
m
S
m
= ϕ
T
m
Bϕ
m
− ϕ
T
m
BϕΣϕ
T
Bϕ
m
Q
m
= ϕ
T
m
By − ϕ
T
m
BϕΣϕ
T
By
http://www.ma-xy.com 30 http://www.ma-xy.com
http://www.ma-xy.com
第一章 回归模型 1.4 RVM 稀疏向量机
Step11. 终止条件，如果算法不终止，则返回 Step2。
基本的稀疏贝叶斯模型至此已经结束了，上面介绍了 (回归问题) 稀疏贝叶斯模型以及模型
的参数估计方法：第二类极大似然估计 (证据近似)、EM 算法以及快速序列算法。关于用 RVM
处理分类 (二分类、多分类) 问题，可以参考 Tipping 的论文，也可以参考 PRML 书籍。下面介
绍 RVM(基) 核函数的改进方法 - 多核方法，此方法在许多核方法中同样适用，比如前面介绍的
SVM。
1.4.5
核方法扩展
-
多核方法
值得一提的是，RVM 的核可以是任意核函数，不像 SVM 中核函数需要满足“Mecer”定理，
这使得我们能够较自由的选择核函数。先来看下面的核函数
ϕ(x) =
K(x
1
, x
1
) K(x
1
, x
2
) . . . K(x
1
, x
M
)
K(x
2
, x
1
) K(x
2
, x
2
) . . . K(x
2
, x
M
)
.
.
.
.
.
.
K(x
N
, x
1
) K(x
N
, x
2
) . . . K(x
N
, x
M
)
N×M
= [ K(·, x
1
), K(·, x
2
), . . . , K(·, x
M
)]
T
≡ [ϕ
1
(x), ϕ
2
(x), . . . , ϕ
M
(x)]
T
其中：
K(x
i
, x
j
) = exp
−
||x
i
− x
j
||
2
2θ
2
可以看出，上面的核函数 K(x
i
, x
j
) 是单一的，不存在其它形式的核函数，下面考虑将多个
核函数混合使用，并且由于 RVM 对核函数没有限制，组合核函数更是多种多样的。常见的高斯
核是经典的局部核，多项式和是经典的全局核。现在将二者结合。定义一个组合核：
K(x
i
, x) = β
1
G
1
(x
i
, x) + β
2
G
2
(x
i
, x) + β
3
P (x
i
, x)
其中：
G
1
(x
i
, x) = −
||x
i
− x||
2
θ
2
1
G
2
(x
i
, x) = −
||x
i
− x||
2
θ
2
2
P (x
i
, x) =
x
i
x
θ
2
3
+ 1
2
θ
1
, θ
2
, θ
3
为核校正后得到的核参数，β
1
, β
2
, β
3
为待求的组合权重。
提到组合核，自然想到组合投资、组合预测等等，其基本工作就是求解组合因子的权重 β
i
。
http://www.ma-xy.com 31 http://www.ma-xy.com
http://www.ma-xy.com
1.4 RVM 稀疏向量机 第一章 回归模型
更一般的，考虑有 M 个核函数进行组合，定义其组合核为
K(x
i
, x) =
M
j=1
β
j
K
j
0 ⩽ β
j
⩽ 1
j
β
j
= 1
其中：K
j
为基本核函数，M 为核函数的个数，β
j
为 K
j
的权重系数。如果在 M 个核函数中，
∀i, j ∈ M，K
i
和 K
j
是相同类型的核，比如：M 个核函数都是高斯核，只是核函数 θ 不同，则
称组合核 K 为同构多核 (同样结构的核)；同样，可以定义异构多核 (前面的示例即为异构 3 核)。
如果给出了 β
j
, K
j
，K 也就确定了，线性函数就变为
y = wK + ε
现在的问题是：应该如何确定组合核 K，或者在给定基本核 K
1
, K
2
, . . . , K
M
后，如何确定组合
权重 β？首先，定义核函数的内积：
定义 (核函数的内积) 给定样本 D = {x
i
, y
i
}
N
i=1
和两个核函数 K
1
, K
2
，K
1
, K
2
所形成的内
积为
K
1
, K
2
F
=
N
i,j=1
K
1
(x
i
, x
j
)K
2
(x
i
, x
j
)
定义 (核校正) 定义核校正为
A(D; K
1
, K
2
) =
K
1
, K
2
F
K
1
, K
1
F
K
2
, K
2
F
由核校正的定义，我们选定的组合核函数 K 和理想核函数 yy
T
之间的核校正值为
A(D; K, yy
T
) =
K, yy
T
F
K, K
F
yy
T
, yy
T
F
=
K, yy
T
F
n
K, K
F
注：这里我们讨论的分类问题 y
i
∈ {−1, +1}。对于分类问题 y(x) 而言，假设分类标签已知，那
么最理想的核函数是 K(x
i
, x
j
) = y(x
i
)y(x
j
)，因为
y(x
i
) = y(x
j
) ⇒ K(x
i
, x
j
) = 1
y(x
i
) = y(x
j
) ⇒ K(x
i
, x
j
) = −1
若 A 的值越高，则表明 K 与 yy
T
越接近 (相似)，从而最佳组合核 K 的确定过程为
max
β
A(β) = A(D; K(β), yy
T
) =
K, yy
T
F
n
K(β), K(β)
F
s.t.
0 ⩽ β
j
⩽ 1
j
β
j
= 1
(1.17)
http://www.ma-xy.com 32 http://www.ma-xy.com
http://www.ma-xy.com
第一章 回归模型 1.4 RVM 稀疏向量机
如果 β
j
→ 0，则证明 β
j
对应的基本核函数 K
j
不适合数据 D。下面来分析上面建立的优化模型
(1.17)
max
β
A(β) =
k
β
k
y
T
K
k
y ≡ A
n
kl
β
k
β
l
K
k
, K
j
F
≡ B
s.t.
0 ⩽ β
j
⩽ 1
j
β
j
= 1
(1.18)
我们要在约束下最大化目标 A(β) = A/B，而最大化 A/B 可以表述为 max A 并且 min B，
于是，将模型 (1.18) 改为
max
β
k
β
k
y
T
K
k
y −
kl
β
k
β
l
K
k
, K
j
F
s.t.
0 ⩽ β
j
⩽ 1
j
β
j
= 1
(1.19)
将上述模型 (1.19) 写为向量形式，且转化为最小化问题，有
min
β
B
T
QB − C
T
B
s.t.
0 ⩽ B ⩽ 1
e
T
B = 1
(1.20)
其中：B = [β
1
, β
2
, . . . , β
M
]
T
，C = [y
T
K
1
y, y
T
K
2
y, . . . , y
T
K
M
y]
T
，
Q =
K
1
, K
1
K
1
, K
2
. . .
K
1
, K
M
K
2
, K
1
K
2
, K
2
. . .
K
2
, K
M
.
.
.
.
.
.
.
.
.
.
.
.
K
M
, K
1
K
M
, K
2
. . .
K
M
, K
M
上述 (1.19) 或者 (1.20) 是一个凸二次规划 (Convex QP) 问题，可以用一般的二次规划算法
求解，可以参考优化部分相应的章节。现在，将这个凸二次规划问题转化为二阶锥规划 SOCP。
令 B
T
QB − C
T
B = t，则凸二次规划变为如下凸二次约束优化问题
min
β
t
s.t.
0 ⩽ B ⩽ 1
e
T
B = 1
B
T
QB − C
T
B ⩽ t
http://www.ma-xy.com 33 http://www.ma-xy.com
http://www.ma-xy.com
1.4 RVM 稀疏向量机 第一章 回归模型
上式是一个以 t, β 为变量的凸二次约束优化问题 (QCQP)，其中：Q 为实对称矩阵 Q = S
T
ΛS =
(
√
ΛS)
T
(
√
ΛS)。于是，上面的 QCQP 可以写为
min
β
t
s.t.
0 ⩽ B ⩽ 1
e
T
B = 1
||DB||
2
⩽ 1 + t + C
T
B
(1.21)
引理 ∀v ∈ R
n
，x, y ∈ R
+
，如果 ||v||
2
⩽ xy，则
2v
x − y
⩽ xy
证明 因为
||v||
2
⩽ xy
4||v||
2
⩽ 4xy
(x − y)
2
+ 4||v||
2
⩽ (x + y )
2
所以
2v
x − y
⩽ xy
□
由上面的引理，可以将核优化的 QCQP 问题 (1.21) 约束中的 ||DB||
2
⩽ 1 + t + C
T
B 改写
min
β
t
s.t.
0 ⩽ B ⩽ 1
e
T
B = 1
2DB
1 − t −C
T
B
⩽ 1 + t + C
T
B
(1.22)
上面的优化模型 (1.22) 即为一个二阶锥规划 SOCP，可以用一般的优化软件进行求解，如：
SeDuMi、YALMIP 和 MATLAB 自带的优化工具箱。
注：二阶锥规划 SOCP 是一个非光滑非线性的凸规划问题，许多优化问题 (如：线性规划 LP、
凸二次规划 QP、凸二次约束规划 QCQP 等) 都可以转化为 SOCP 问题，因此，SOCP 是最优
化学科的一个研究热点。SOCP 一般采用内点法进行求解，常用的工具箱有 SeDuMi 和 Yalmip，
SeDuMi 求解 SOCP 的迭代次数为 O(
√
M log(1/ε))，其中：ε 为求解误差。关于 SOCP 更详细
的理论，可以参考本书最优化部分。
http://www.ma-xy.com 34 http://www.ma-xy.com
http://www.ma-xy.com
第一章 回归模型 1.5 MATLAB 回归简介
1.5 MATLAB 回归简介
1.5.1 MATLAB 回归命令
下面，将介绍 MATLAB 的可用于回归问题的一些命令。
(1) regress
可用于多元线性回归，其命令格式为
[b,bint,r,rint,stats] = regress(y,x,alpha)
其中：b 为回归系数，即权重；bint 为回归系数的 (1-alpha)×100% 置信区间；r 为残差；rint 为
残差的 (1-alpha)×100% 置信区间；stats 为结构体，包含判决系数 R
2
、样本的 F 统计量值、检
验的 p 值以及误差方差 σ
2
的估计值。值得一提的是，如果我们要估计模型中的常数项，则 x 的
第一列设置为全 1 列。
(2) regstats 可用于多重线性或广义线性回归，其命令格式为
[stats] = regstats(y,x,model,whichststs)
其中：stats 为结构体变量，返回所有诊断统计量 (关于线性回归模型的诊断，可以参考一般的计
量经济学课本)，如果省略输出 stats，会有交互式 GUI 界面；x 会自动增加全 1 列，即模型中默
认存在常数项；model 可选的参数值有：’linear’(有常数项)、’interaction’(模型中有常数项、线性
项和交叉项)、’quadratic’(常数项、线性项、交叉项和平方项)、’purequadratic’(常数项、线性项
和平方项)；whichstats 用于指定输出的统计量，示例’leverage’,’standres’。
(3) robustt 用于稳健回归，应该是用到加权最小二乘估计法，模型中参数受异常值的影响
较小。其命令格式为
[b,stats] = robustt(x,y,wfun,tune,const)
其中：b 为回归系数；stats 为结构体，包含 ols.s(ˆσ by RMSE)、robust.s(ˆσby robust)、mod.s(对
残差绝对值中位数计算 σ)、s(ˆσ 最终值)、se(模型系数的标准误差)、t(样本的 t 统计量值，是 b
和 se 的比值)、p(t 检验的 p 值)、corb(b 的协方差矩阵估计值)、coeurr(b 的相关系数估计值)、
w(robust 权重向量)、h(ols 的中心化抽样值向量)、dfe(误差的自由度)、R(x 的 QR 分解中的 R
因子)；x 为 n × p 数据；y 为 n × 1 数据；wfun 为加权函数可以列表形式展示一下；tune 为常
数 r =
resid
tune×s×
√
1−h
；const 是否添加常数项，’on’ 表示添加。
(4) ridge 就是 ridge 回归，其命令格式为
[b] = ridge(y,x,k,scaled)
其中：y 是 n ×1 向量；x 是 n ×p 矩阵；k 为参数向量，如果可有 m 个元素，则 b 是一个 p ×m
元素，默认情况下是均值为 0 方差为 1 的正态分布，模型不包含常数项；scaled 用于指定是否包
含常数项，scaled=0 表示包含常数项。
(5) lasso 用于 lasso 回归，其命令格式为
[B,FitInfo] = lasso(x,y,params)
其中：B 是模型系数，是一个 p ×l 矩阵，l 是 lambda 的数量；params 可写的参数和参数值有：
’Alpha’, 取值为 [0, 1]，当取值为 0 时表示 ridge 回归，取值为 1 时表示 lasso 回归，默认为 1；
’CV’ 用于 k 折交叉验证等等。其示例为
1 load co rbig
2 x = [ Ac c ele r ati o n Displacement Hovsepouer Weight ] ;
http://www.ma-xy.com 35 http://www.ma-xy.com
http://www.ma-xy.com
1.6 RVM 工具箱 - SB2 第一章 回归模型
3 [ b , f i n i n f o ] = l a s s o ( x ,MPG, ’CV’ ,1 0 ) ;
4 l a s so Pl o t (b , ’ f i t i n f o ’ , ’ plotType ’ , ’ lambda ’ , ’ XScale ’ , ’ lo g ’ ) ;
5 %构 建 E l a s t i c Net
6 nonam = ~any ( is nan ( [ x ,MPG] ) ,2 ) ;
7 Xnoneam = x(nonam , : ) ;
8 MPGnonam = MPG(nonam , : ) ;
9 [ ba , f i t i n f o a ] = l a s s o (Xnonam,MPGnonam, ’CV’ ,10 , ’ alpha ’ ,0 . 5 ) ;
10 pnames = { ’ Acc e l e ra io n ’ , ’ Displacement ’ , ’ Horsepower ’ , ’Weight ’ };
11 l a s so Pl o t ( ba , f i t i n f o a , ’ PlotType ’ , ’ lambda ’ , ’ XScale ’ , ’ l o g ’ , ’ PredretorNames ’ , ’ pnames ’ ) ;
12
(6) PLSregress 用于偏最小二乘回归，其命令格式为
[XL,YL,Xs,Ys,Beta,PctVar,Mse,Stats] = PLSregress(x,y,ncomp,params)
其中：x 为 n×p 矩阵；y 为 n×1 向量；ncomp 为 Number of PLS components；XL 为 p ×ncomp
矩阵；YL 为 m × ncomp 矩阵；params 用来设置参数和参数值，可用参数有：’CV’、’mcreps’
和’options’。另外，MATLAB 有用于 PLS 的工具箱 PLStool。
(7) nlint 用于非线性回归，利用 Levenberg-Marquardt 算法求解最小二乘优化问题，其命
令格式为
[beta,r,J,covb,mse] = nlint(x,y,fun,beta0,options)
其中：J 为雅可比矩阵；mse 为误差方差的均方误 mse；fun 用于指定非线性函数，如 @(x)y = ax
2
；
beta0 为模型参数的初始值；options 用于指定参数和参数值，如’Display’,’o’,’MaxIter’,1000,
’TolFun’(残差平方和-模型的终止容限),1e8,’TolX’(参数估计值的终止容限),1e8,’FunValCheck’,’on’-
检验 y 的无效值 NaN 或者 inf,’Robust’,’on’-用加权最小二乘稳健估计,’WgtFun’-Weightfun,’Tune’-
稳健拟合的调节参数。另外，在非线性拟合结束之后，可以用 nlpredci 求参数的置信区间。
1.5.2 MATLAB 回归示例
1.6 RVM 工具箱 - SB2
todo: 待补充。。。
1.7 非参数回归
上面介绍的回归方法都是参数回归，例如线性回归、广义线性回归以及贝叶斯回归。下面，
我们将介绍一些常用的非参数回归 (拟合) 方法，包括：基于局部拟合的局部常系数核权重拟合、
局部多项式核权重拟合和 k 邻近核权重估计；基于样条方法的光滑样条拟合、多项式样条拟合和
罚样条；基于小波方法的小波回归以及正交序列回归。之所以说回归拟合，是因为其本质工作是
一样的，并且，由于密度估计和回归估计是相似的，所以下面介绍的这些非参数方法可以推广到
非参数密度估计。
http://www.ma-xy.com 36 http://www.ma-xy.com
http://www.ma-xy.com
第一章 回归模型 1.7 非参数回归
1.7.1 核估计
局部常系数核权重拟合
设 X 为输入变量，Y 为输出变量，(x
i
, y
i
)
n
i=1
为 n 个样本，并且为了简便，我们设 X, Y ∈ R。
现在考虑 X 到 Y 的回归问题，即求 Y = f(X), f : R → R。在前面的参数回归中，我们是在正
态假设的前提下，用 ML 等方法求解 y = f (x)，这是一种设定 f 的形式，在正态假设下求 f 中
参数 θ 的方法。在接下来的非参数回归方法中，我们扔掉正态前提。回归 (拟合) 问题，从局部
来看，只要估计每个 x
i
对应的 ˆy
i
= f(x
i
) 即可，这类似于一个计分方法，比如：假设有如下表
(1.3) 的数据
表 1.3: 非参数回归模拟数据
样本 x y
1 0.5 1
2 0.4 1.3
3 0.51 1.1
4 0.48 0.9
如果要估计 x
1
处的 f (x
1
)，一种合理的想法是：
w
i
y
i
，其中：w
i
是权重，取决于 x
i
与 x
1
的距离，w
i
距离 x
1
近则 w
i
大。我们可以设置权重 w
i
为密度函数 (例如正态分布密度)，于是
有
ˆ
f(x) =
n
i=1
w
i
(x)y
i
。w
i
是 x 的函数，与 x 的位置有关，并且所有样本点皆参加 x 处样本
值的估计。统计学家用核作为权重，称为非参数核回归，其形式如下
w
i
(x) =
K
x
i
−x
h
n
n
i=1
K
x
i
−x
h
n
∈ [0, 1]
其中：h
n
为窗宽，与样本量 n 有关。常见的核函数有：
1. 均匀核：K(u) =
1
2
I(|u| ⩽ 1)；
2. 高斯核：K(u) =
1
√
2π
e
−
u
2
2
；
3. Epanechnikov
核：
K
(
u
) =
3
4
(1
−
u
2
)I(|u| ⩽ 1).
在 K(u) 为均匀核时，
ˆ
f
K
(x) 就是落在邻域 [x − h
n
, x + h
n
] 中的样本 x
i
对应的 Y
i
的算术平均
值。这是 Nadaraya 和 Watson 于 1964 年分别给出的方法，我们记此回归为
ˆ
f
NW
(x)，∀x
i
，只
要给定 K(·) 和 h
n
，即可求出其估计值
ˆ
f
NW
(x
i
)。
对于一般的回归问题，我们要讨论回归估计得好坏，以及估计量的统计性质 (相合性、强相
合性、渐进偏差、渐进正态等)。如果模型中有外来参数 K(·), h
n
，我们要讨论如何选择它们 (如
何设置准则进行比较)。
下面，先来讨论模型拟合结果的好坏。制定一个准则，优良的 K(·), h
n
应该在一定的准则下
使模型拟合效果最好。这里我们仅讨论 h：
http://www.ma-xy.com 37 http://www.ma-xy.com
http://www.ma-xy.com
1.7 非参数回归 第一章 回归模型
定义风险 (均方误差) 为
R(h) = E
1
n
n
i=1
(
ˆ
f(x
i
) − f (x
i
))
2
风险 R(h) 越小 h 越好。但是，一般情况下 f(x) 是未知的，于是，我们用观测值 y
i
来代替 f(x
i
)。
但是这样做会低估风险，因为我们在估计
ˆ
f(x
i
) 时用了一次 y
i
，之后再用一次的话，显然低估了
风险。我们采用一种交叉验证的方法来选取 h
h
∗
= min
h
CV (h) =
1
n
n
i=1
[Y
i
−
ˆ
f
−i
(x
i
)]
2
=
1
n
2
h
n
i=1
n
j=1
K
x
i
− x
j
h
−
2
n
(
n
−
1)
h
n
i=1
n
j=i
K
x
i
− x
j
h
其中：
ˆ
f
−i
(
x
i
)
表示在估计
x
i
处的
ˆ
f
(
x
i
)
时，样本值
y
i
不加入计算。这里的
h
∗
是可行的，但其
最优化的推导求解并不是容易的。此外，对于 h 的选取，还有经验法、插值法、AIC 方法以及
GCV 方法等。
我们不加证明的给出估计量
ˆ
f
NW
(x) 的性质：
1. 逐点矩相合性：
ˆ
f(x)
P
−→ f(x);
2. 全局相合性：E(|
ˆ
f(x) − f(x)|
p
) → 0;
3. 几乎处处强相合性：
ˆ
f(x) → f (x), a.s;
4. 一致强相合性：sup |
ˆ
f(x) − f(x)| → 0, a.s;
5. 渐进正态性：
√
nh
n
[
ˆ
f(x) − E(
ˆ
f(x))]
·
∼ N(0, r(x));
6.
√
nh
n
[
ˆ
f(x) − f(x)]
D
−→ N (0, r(x));
上述性质的详细证明可以参考《现代非参数统计》薛留根 P171 和《现代非参数统计》L. 沃
塞曼 P61。此外，估计量还具有如下渐进无偏
lim
n→∞
E(
ˆ
f(x
i
)) = f(x
i
)
其渐进收敛速度为
E(
ˆ
f(x
i
)) − f (x
i
) = c
1
h
2
+ o(h
2
) + O((nh)
−1
)
其中：
c
1
=
1
2
f
′′
(x
i
) +
ˆ
f
′
(x
i
)f
′
(x
i
)
f
′
(x
i
)
t
2
K(t)dt
NW 估计在边界点 x(1), x(n) 处的估计偏差较大，为此，1984 年 Gasser-Niuller 提出一种新
的核估计
ˆ
f
GM
(x) =
n
i=1
s
i
s
i−1
K
x
i
− x
h
n
dt
y
i
其中：s
i
=
x(i)−x(i+1)
2
，x(0) = −∞，x(n + 1) = ∞。
http://www.ma-xy.com 38 http://www.ma-xy.com
http://www.ma-xy.com
第一章 回归模型 1.7 非参数回归
局部多项式核权重拟合
不难看出，上面我们假设了 f (x) 在任意点 x
0
处的邻域内是一个常数，核回归是通过加权最
小二乘方法得到 x
0
处的估计值。设 x
0
出的常数为 a =
ˆ
f(x
0
)
min
a
n
i=1
(y
i
− a(x
0
))
2
w
i
(x
0
) (1.23)
但实际上，函数 f(x) 在 x
0
的邻域内并非一个常数，由 Tayler 公式有
f(x) ≈ f (x
0
) + f
′
(x
0
)(x − x
0
) +
f
′′
(x
0
)
2!
(x − x
0
)
2
+ ··· +
f
p
(x
0
)
p!
(x − x
0
)
p
=
p
j=0
β
j
(x
0
)(x − x
0
)
j
可以看出，函数 f(x) 在 x
0
的邻域内不止可以设置为常数，而且可以设置为 p 次多项式，于是
我们改进目标 (1.23)，有
min
β
n
i=1
y
i
−
p
j=0
β
j
(x
0
)(x
i
− x
0
)
j
2
w
i
(x
0
) (1.24)
上式即所谓的局部多项式核权重拟合。其中：参数为 p, h
n
, K(·)，β
j
(x
0
) 为所求，并且通常
情况下，我们选取 p 为奇数。在 p, h
n
, K(·) 确定之后，每给一个点 x
i
，可以通过加权最小二乘
方法得到 β
j
(x
0
)，于是会有 x
0
处的估计值
ˆ
f(x
0
)。为了求 β(x
0
)，记
X
x
0
=
1 x
1
− x
0
. . .
(x
1
−x
0
)
p
p!
1 x
2
− x
0
. . .
(x
1
−x
0
)
p
p!
.
.
.
.
.
.
.
.
.
.
.
.
1 x
n
− x
0
. . .
(x
n
−x
0
)
p
p!
令 W
x
0
为 n × n 对称矩阵，对角为 diag(w
i
(x
0
))，则目标 (1.24) 可改写为
(y − X
x
0
β)
T
(y − X
x
0
β)
于是有
ˆ
β(x
0
) = (X
T
x
0
W
x
0
X
x
0
)
−1
X
T
x
0
W
x
0
y
上面介绍了两种核思路：一种是局部常数，一种是局部多项式，那么，对于不同的核估计，应该
如何做比较呢？此类问题是统计学中最小最大风险，可以参考相关的数理统计书籍或者《现代非
参数统计》P227。
k 邻近核权重估计
在前面介绍的两种核估计方法中，对点 x
0
处的函数值进行估计时，所有样本点 (x
i
, y
i
)
n
i=1
都参与计算，这样的方法计算量大而且这种思路也不是很合理，其实用 x
0
附近的几个点来计算
http://www.ma-xy.com 39 http://www.ma-xy.com
http://www.ma-xy.com
1.7 非参数回归 第一章 回归模型
也就足够了。为此，引出 k 邻近核权重估计，我们用 x
0
附近的 k 个样本点来进行估计，有
ˆ
f
NN
(x
0
) =
k
i=1
w
i
(x
0
)y
i
其中：权重为
w
i
(x
0
) =
K
x
i
−x
0
max |x
i
−x
0
|
k
i=1
K
x
i
−x
0
max |x
i
−x
0
|
注：k 个样本点可以是变动的 k
x
，这样有利于处理一些离群值，并有助于控制计算量。同局部常
数核估计一样，只不过这里选取了 k 个样本点，窗宽为 max |x
i
−x
0
|。其改进方向有：k 邻近多
项式拟合以及权重 w
i
(x) 形式的改进。
我们不加证明的给出 k 邻近核权重估计的性质：
(1) 渐进无偏性。
ˆ
f
NN
(x) 的渐进偏差为
u(k)
8f(x)
3
[
ˆ
f
′′
f + 2
ˆ
ff
′
(x)]
k
n
2
渐进方差为
2
σ
2
(x)
k
R(k)
(2) 相合性。
ˆ
f
NN
(x) → f(x), a.s
(3) 渐进正态性。
nh
n
[
ˆ
f
NN
(x) − f (x)]
K
−→ N(0, r
2
(x))
其中：r
2
(x) = V ar(y|x)
∞
−∞
K
2
(u)du。
1.7.2 样条拟合
光滑样条拟合
所谓样条就是拟合分段多项式，用样条逼近函数曲线的方法称为样条方法，逼近函数的导数
的间断点称为“节”。样条实际上可用来逼近任何光滑函数，至少当节点数据充分大时，可以做到
这一点。为引出样条方法，我们先回到最本质的问题：有待求函数 f(x) 和样本数据 (x
i
, y
i
)，我
们的目标是
ˆ
f
∗
= min
ˆ
f
n
i=1
(y
i
−
ˆ
f(x
i
))
2
http://www.ma-xy.com 40 http://www.ma-xy.com
http://www.ma-xy.com
第一章 回归模型 1.7 非参数回归
上述问题是一个泛函问题，找一个最佳估计函数 f
∗
，使其与 y 的离差平方和最小，如果不加任
何限制，即
ˆ
f 可以是任何函数，我们可以找到无数个最优解
ˆ
f
∗
。但这样做对数据是过拟合的，一
点误差都没有是没有意义的。我们可以想象这些最优解 (内插函数){
ˆ
f∗} 的特点：很多是非连续
的，不可导的，它们恰好在 x
i
处的函数值为 y
i
。而我们需要的 f (x) 往往要具有一定的可导可微
性，于是，我们使用下面这种做法
min
ˆ
f
M(λ) =
n
i=1
(y
i
−
ˆ
f(x
i
))
2
+ λJ(
ˆ
f(x))
通常取 J
k
(
ˆ
f(x)) =
1
0
f
(k)
(x)dx，k 在很多情况下为 2，反应了
ˆ
f 的光滑程度。λ 越大函数
ˆ
f 越
光滑，当 λ → 0 时，
ˆ
f 收敛到最小二乘直线。这种在优化对象后面加一个惩罚的做法称为正则
化，我们前面介绍过 ridge 和 lasso，这样得到的估计量记为
ˆ
f
ss
(x)。在样条方法中，这种做法称
为光滑样条估计，当 0 < λ < ∞ 时，
ˆ
f 在内插函数和一次函数之间波动，一般的解为三次样条。
三次样条在每个点 x
i
有节，在 [x(i), x(i + 1)] 为三次多项式，且 2 阶可导，并且，三次样条在
[x(i), x(i + 1)] 是唯一解。例如：我们设定三次多项式为
B
3
(x) = a
0
+ a
1
x
1
+ a
2
x
2
+ a
3
x
3
并设定共有 5 个样本点 x
i
, y
i
，样本的顺序统计量为 x(i)，我们共有 4 段 [x(i), x(i+1)], i = 1, 2, 3, 4，
每一段上都是一个三次多项式，并且该三次多项式满足这一段的端点条件，我们共有 4 个 3 次
多项式，其示意图如图 (1.5) 所示
图 1.5: 三次样条拟合示意图
我们用 a
ij
表示第 i 个三次多项式的 x
j−1
项的系数，其待求系数如表
表 1.4: 三次样条拟合系数表
项 常数项 x
1
x
2
x
3
1 a
11
a
12
a
13
a
14
2 a
21
a
22
a
23
a
24
3 a
31
a
32
a
33
a
34
4 a
41
a
42
a
43
a
44
上面的拟合虽然在节点处很完美，但是当样本量较少时，则为完全拟合且函数整体波动较大，
当样本量比较大时，整体波动较小，但计算量相对较高。上面的三次样条设置为 B
3
(x)，并且有
4 个三次样条，称这 4 个三次样条为样条基函数。现在设 B 样条为 B(x)，拟合函数为
ˆ
f(x) =
N
j=1
ˆ
β
j
B
j
(x)
http://www.ma-xy.com 41 http://www.ma-xy.com
http://www.ma-xy.com
1.7 非参数回归 第一章 回归模型
这样，只需要找到系数
ˆ
β
j
即可。将目标函数写为矩阵形式，有
(y − Bβ)
T
(y − Bβ) + λβ
T
Ωβ
其中：B
ij
= B
j
(x
i
)，Ω
jk
=
1
0
B
′′
j
(x)B
′′
k
(x)dx。由矩阵运算，有
ˆ
β = (B
T
B + λΩ)
−1
B
T
y
记 L = B(B
T
B + λΩ)
−1
B
T
为帽子矩阵，则
ˆ
f(x) = Ly
多项式样条拟合
由上面的启发，我们直接用样条函数 B(x) 来逼近 f (x)，就像以前直接用 y = ax
2
来逼近
f(x) 一样，设置目标为
min
β
n
i=1
[y
i
− β
T
B(x
i
)]
2
其中：B(x) = (B
1
(x), B
2
(x), . . . , B
q
(x))
T
是一个样条基。现在，问题来了，不同节点位置，不同
节点数量会导致拟合效果截然不同，如此，我们应该如何确定节点数目 J 和节点位置 t
j
？
罚样条
惩罚样条是稳健估计，其光滑参数可以用交叉验证方法来求解。下面，我们来介绍两种罚样
条：幂基样样条和 B 样条。
(1) 幂基样条。幂基样条的优化目标为
min
β
n
i=1
[y
i
− β
T
B(x
i
)]
2
+ λ
J
j=1
β
2
q+j
其中：β 为所有样条系数，B 为幂基样条基，J 为节点数，q 为样条多项式阶数
1, x, . . . , x
q
, (x − t
1
)
q
+
, . . . , (x −t
J
)
q
+
这里的 a
+
= max{0, a}。将上面的优化模型写为矩阵形式，有解
ˆ
β = (B
T
B + λD
q
)
−1
B
T
y
其中：B = (B
T
(x
1
), . . . , B
T
(x
n
))
T
，y = (y
1
, y
2
, . . . , y
n
)
T
，D
q
= diag(0
q +1
, 1
j
)，
ˆ
f = B
ˆ
β = B(B
T
B + λD
q
)
−1
B
T
y
(2)B 样条。B 样条的优化目标为
min
β
n
i=1
[y
i
− β
T
B(x
i
)]
2
+ λ
q+J
j=k+1
(∆
k
β
j
)
2
http://www.ma-xy.com 42 http://www.ma-xy.com
http://www.ma-xy.com
第一章 回归模型 1.7 非参数回归
其中：∆ 为微分算子，
∆
k
β
j
=
k
l=0
(−1)
l
C
l
k
β
j−l
节点为等距节点。我们记 B = (B
T
(x
1
), . . . , B
T
(x
n
))
T
，y = (y
1
, y
2
, . . . , y
n
)
T
，D
k
是 (q + J −k) ×
(q + J) 矩阵
D
k
=
(−1)
0
C
0
k
. . . (−1)
k
C
k
k
0 . . . 0
0 (−1)
0
. . . (−1)
k
C
k
k
. . . 0
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
0 . . . 0 (−1)
0
C
0
k
. . . (−1)
k
C
k
k
其中：C
r
k
=
k
r
=
k!
r!(k−r)!
。于是目标可写为矩阵形式，有
(y − Bβ)
T
(y − Bβ) + λβ
T
D
T
k
D
k
β
得到参数估计为
ˆ
β = (B
T
B + λD
T
k
D
k
)
−1
B
T
y
且
ˆ
f(x) = B
ˆ
β
关于光滑参数 λ 的选取，可以用 CV 和 GCV 进行确定。其去 1 交叉验证为
CV (λ) =
1
n
n
i=1
y
i
−
ˆ
f
λ
(x
i
)
1 − h
ii
2
其中：h
ii
是帽子矩阵 L(λ) 的对角线元素，幂基的帽子矩阵为 B(B
T
B + λD
q
)
−1
B
T
，B 基的帽
子矩阵为 B(B
T
B + λD
T
k
D
k
)
−1
B
T
。去 1 广义交叉验证为
GCV (λ) =
1
n
n
i=1
y
i
−
ˆ
f
λ
(x
i
)
1 − n
−1
tr(L(λ))
其中；tr(L(λ)) 为 L(λ) 的迹。
1.7.3 正交级数回归
先来看正交级数回归。定义正交基 {φ
i
(x)}
1
0
φ
i
(x)φ
j
(x)dx =
1 i = j
0 i = j
http://www.ma-xy.com 43 http://www.ma-xy.com
http://www.ma-xy.com
1.7 非参数回归 第一章 回归模型
由于待求函数 f(x) 可以用 {φ
i
(x)} 的 Fourior 级数表示
f(x) =
∞
j=1
a
j
φ
j
(x)
于是对 f(x) 的估计就变为对系数 a
j
的估计。我们知道 a
j
为
a
j
=
1
0
f(x)φ
j
(x)dx
对 a
j
最直接的估计为
ˆa
j
=
1
0
y
i
φ
j
(x)dx =
n
i=1
(x
i
− x
i−1
)y
i
φ
j
(x
i
)
在前面的问题说明 (1.1) 当中，我们提到过回归数据有两种形式，一种是实验数据，一种是
统计数据，在这里，我们称之为固定模式和随机模式。
(1) 对固定模式有 0 = x
1
⩽ x
2
⩽ ··· ⩽ x
n
= 1 。当 x
i
为 [0, 1] 上的等距点时，有 a
j
的估计
ˆa
j
=
n
i=1
(x
i
− x
i−1
)y
i
− φ
j
(x
i
)
=
1
n
n
i=1
y
i
φ
j
(x
i
)
以及 f(x) 的估计
ˆ
f(x) =
N
j=1
ˆa
j
φ
j
(x
i
)
其中：N 为基组中基的个数。
(2) 对随机模式，样本 x
i
, y
i
独立同分布，
f(x) ≈
k
j=1
a
j
φ
j
(x)
a
j
=
1
n
n
i=1
f(x
i
)φ
j
(x
i
)
其估计量为
ˆa
j
=
1
n
n
i=1
y
i
φ
j
(x
i
)
ˆ
f(x) =
k
j=1
η(ˆa
j
)φ
j
(x)
其中：η(·) 为门限函数，当 |a| ⩾ δ
n
时，η(a) = a，否则为 0。关于基个数 k 以及门限 δ
n
的选
取，这里不做介绍。
http://www.ma-xy.com 44 http://www.ma-xy.com
http://www.ma-xy.com
第一章 回归模型 1.7 非参数回归
关于正交基，在前面的偏微分方程章节的谱分析中介绍了一些，比如 Legendre 正交多项式和
Chebyshev 正交多项式等等，这里不再介绍。值得一提的是，这些正交多项式都是在一定的自变
量 x 范围内正交的，比如 Legendre 在 [−1, 1] 上正交。对于实验数据和统计数据，如果 x ∈ [a, b]，
我们可以做变换 z =
2x−a−b
b−a
∈ [−1, 1] 来把区域 [a, b] 投射到 [−1, 1]。
1.7.4 小波核估计
局部常数小波核权重拟合
最流行的正交基 {φ
i
(x)} 是小波正交基，φ
j
是由父小波和母波的转换构成的，关于小波的
内容，可以参考《统计建模的小波方法》Brani Vidakovic 著田铮译。
(1) 对固定模式。设 x
i
非随机，且 0 = x
1
⩽ x
2
⩽ ··· ⩽ x
n
= 1 ，有
ˆ
f(x) =
n
i=1
y
i
s
i
s
i−1
K
j
(x, y)dy
其中：s
0
= 0 ，s
n
= 1 ，s
i
=
x
i
+x
i+1
2
。一般采用小波 GasserMuller 核估计。
(2) 对随机模式。(x
i
, y
i
) 独立同分布，条件期望 f (x
i
) = E(y
i
|x
i
)，对此，Antoniadis, gregoire 和
Mckegve 建议用 NW 小波形式
ˆ
f(x) =
n
i=1
y
i
K
j
(x, x
i
)
n
i=1
K
j
(x, x
i
)
其中：K
j
(x, y) = 2
j
k∈Z
ϕ(2
j
x −k)ϕ(2
j
y −k)，ϕ 为尺度函数，j 相当于核估计中的窗宽，可以用
CV 和 GCV 等方法选取。
局部多项式小波核权重拟合
在核估计中的局部多项式核权重估计中，我们将核改为小波核，有
min
β
n
i=1
y
i
−
p
j=0
β
j
(x
i
− x)
K
j
(x
i
, x)
关于非参数回归方法，我们先介绍到这里。还有一些其它的估计方法，比如：¬局部自适应
回归样条估计；趋势滤波的自适应分段多项式估计。对于趋势滤波的自适应分段多项式估计，
2009.Kim 对线性趋势滤波 (k = 1) 采用 primal - dual 内点法求解，网址
®
有相应的 MATLAB 和
C 代码。2011.Tibshirani 和 Taylor 采用轨道算法对其进行求解，R 包 genlasso 中的 trendlter
函数是可用的，在 CRAN 库中可以找到。R 的非参数回归命令 ksmooth 做 NW 核回归光滑；
loess 用于局部多项式回归拟合；lowess 用于局部加权描点光滑。
®
http://stanford.edu/boyd/11 _tf
http://www.ma-xy.com 45 http://www.ma-xy.com
http://www.ma-xy.com
参考文献
[1] Duchi. Adaptive subgredient methods for online learning and stochiastic optimization. 2011.
[2] Nesterov. A method for unconstrained convex minimization problem with the rate of con-
vergence o(1/k). 1983.
[3] Qian. On the momentum termin gradient descent learning algorithms. 1999.
[4] Tipping. Sparse bayesian learning and the relevance vector machine. 2001.
[5] Tipping. Fast marginal likelihood maximination for sparse bayesian model. 2003.
[6] Zeiler. Adadelta:an adaptive learning rate method. 2012.
47
