一、牛顿方法
梯度下降数学上较清晰地给出了迭代参数的方法。但有一种通常情况下比它更快的算法,称为牛顿方法。牛顿方法的收敛速度理论上为二次收敛,即迭代前误差为0.1量级,迭代一次后误差就为0.01量级,再迭代一次0.0001,再一次0.0000001,即迭代速度非常快。但是,这仅仅理论上接近正确值时成立,实际上会因为一些复杂的数学原因没有这么快,但也比梯度下降快。牛顿方法适用于Logistic回归以及各类广义线性模型,如logistic一般情况下往往数十次就可以收敛
二、一元函数的牛顿方法
有一个关于
θ
\theta
θ的函数
f
(
θ
)
f(\theta)
f(θ),若我们需要求
f
(
θ
)
=
0
f(\theta)=0
f(θ)=0的
θ
\theta
θ,(
θ
\theta
θ)是一维参数,则可以随机选择一个初始点
θ
k
\theta_k
θk,然后计算在
θ
k
\theta_k
θk处函数
f
(
θ
)
f(\theta)
f(θ)的切线
l
=
f
(
θ
k
)
+
f
(
θ
k
)
′
(
θ
−
θ
k
)
l=f(\theta_k)+f(\theta_k) ^{'}(\theta-\theta_k)
l=f(θk)+f(θk)′(θ−θk),然后令
f
(
θ
k
)
=
0
f(\theta_k)=0
f(θk)=0可以求得
θ
k
+
1
\theta_{k+1}
θk+1。
θ
k
+
1
\theta_{k+1}
θk+1与
θ
k
\theta_k
θk的差值大小由导数定义可得为
f
(
θ
k
)
f
(
θ
k
)
′
\frac{f(\theta_k)}{f(\theta_k)^{'}}
f(θk)′f(θk),即一次迭代后,
θ
k
+
1
\theta_{k+1}
θk+1=
θ
k
\theta_k
θk-
f
(
θ
k
)
f
(
θ
k
)
′
\frac{f(\theta_k)}{f(\theta_k)^{'}}
f(θk)′f(θk),由此可以不断逼近使
f
(
θ
)
=
0
f(\theta)=0
f(θ)=0的
θ
\theta
θ。
在logistic回归中,参数
θ
\theta
θ的选择根据最大似然性原则选择使
l
(
θ
)
=
l
n
(
L
(
θ
)
)
l(\theta)=ln(L(\theta))
l(θ)=ln(L(θ))最大的
θ
\theta
θ。最值处函数导数为0,因此即寻找使
l
(
θ
)
=
l
n
(
L
(
θ
)
)
′
=
0
l(\theta)=ln(L(\theta))^{'}=0
l(θ)=ln(L(θ))′=0的参数
θ
\theta
θ。由上述参数
θ
\theta
θ迭代公式可知,在一元logistic回归中参数迭代为
θ
k
+
1
\theta_{k+1}
θk+1=
θ
k
\theta_k
θk-
l
(
θ
k
′
)
l
(
θ
k
)
′
′
\frac{l(\theta_k^{'})}{l(\theta_k)^{''}}
l(θk)′′l(θk′)
三、多元函数的牛顿方法
微积分中对多元函数求极值时,往往用到Hessian矩阵。多元函数
f
(
θ
1
,
θ
2
,
θ
3
.
.
.
θ
n
)
的
f(\theta_1,\theta_2,\theta_3...\theta_n)的
f(θ1,θ2,θ3...θn)的Hessian矩阵为
[
∂
2
f
∂
θ
1
2
∂
2
f
∂
θ
1
∂
θ
2
.
.
.
∂
2
f
∂
θ
1
∂
θ
n
∂
2
f
∂
θ
2
∂
θ
1
∂
2
f
∂
θ
2
∂
θ
2
.
.
.
∂
2
f
∂
θ
2
∂
θ
n
.
.
.
∂
2
f
∂
θ
n
∂
θ
1
∂
2
f
∂
θ
n
∂
θ
2
.
.
.
∂
2
f
∂
θ
n
2
]
\left[ \begin{matrix}\frac{\partial ^2f}{\partial \theta_1 ^2} & \frac{\partial ^2f}{\partial \theta_1 \partial \theta_2 } & ... &\frac{\partial ^2f}{\partial \theta_1 \partial \theta_n }\\\frac{\partial ^2f}{\partial \theta_2 \partial \theta_1 } & \frac{\partial ^2f}{\partial \theta_2 \partial \theta_2 } &...&\frac{\partial ^2f}{\partial \theta_2 \partial \theta_n }\\...\\\frac{\partial ^2f}{\partial \theta_n \partial \theta_1 } & \frac{\partial ^2f}{\partial \theta_n \partial \theta_2 } & ... & \frac{\partial ^2f}{\partial \theta_n^2}\end{matrix}\right]
⎣⎢⎢⎢⎢⎡∂θ12∂2f∂θ2∂θ1∂2f...∂θn∂θ1∂2f∂θ1∂θ2∂2f∂θ2∂θ2∂2f∂θn∂θ2∂2f.........∂θ1∂θn∂2f∂θ2∂θn∂2f∂θn2∂2f⎦⎥⎥⎥⎥⎤
当Hessian矩阵的行列式
d
e
t
(
H
)
=
0
det(H)=0
det(H)=0时的参数
θ
^
\hat\theta
θ^=
[
θ
1
,
θ
2
.
.
.
θ
n
]
T
[\theta_1,\theta_2...\theta_n]^T
[θ1,θ2...θn]T处取得极值。当Hessian矩阵为负定矩阵(顺序主子式均小于0或特征值均小于0)时,在该参数处取得最大值。
多元函数的牛顿方法迭代公式为
θ
^
k
+
1
\hat\theta_{k+1}
θ^k+1=
θ
^
k
\hat\theta_k
θ^k -
▽
f
H
\frac{\triangledown f}{H}
H▽f=
θ
^
k
\hat\theta_k
θ^k - H-1
▽
f
\triangledown f
▽f
=
θ
^
k
\hat\theta_k
θ^k - H-1 *
[
∂
f
θ
1
,
∂
f
θ
2
,
.
.
.
,
∂
f
θ
n
,
]
T
[\frac{\partial f}{\theta_1},\frac{\partial f}{\theta_2},...,\frac{\partial f}{\theta_n},]^T
[θ1∂f,θ2∂f,...,θn∂f,]T
由于每次迭代都需要计算Hessian矩阵的逆,对于大量特征(几千个特征)时不适用牛顿方法。
四、广义线性模型
1、广义线性模型
在 https://blog.csdn.net/ShiZhanfei/article/details/84345278与 https://blog.csdn.net/ShiZhanfei/article/details/84704299 中,探讨了线性回归与logistic回归。这两种数学模型事实上都可以归为一类更普遍的模型,即广义线性模型。
最小二乘线性回归基于高斯分布的概率形式,即
P
P
P{
y
;
μ
,
σ
y;\mu,\sigma
y;μ,σ}=
1
2
π
σ
e
−
(
y
−
μ
)
2
2
σ
2
\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(y-\mu)^2}{2\sigma^2}}
2πσ1e−2σ2(y−μ)2
logistic回归基于伯努利分布的概率形式,即
P
P
P{y=1;
θ
\theta
θ}=
θ
\theta
θ
广义线性模型基于指数分布族的概率形式,即
P
P
P{y;
η
\eta
η}=
b
(
y
)
∗
e
η
T
T
(
y
)
−
a
(
η
)
b(y)*e^{\eta^T T(y)-a(\eta)}
b(y)∗eηTT(y)−a(η),其中,参数
η
\eta
η称为分布的自然参数,为
θ
T
x
\theta^T x
θTx。T(y)称为充分统计量(统计量是指与样本相关不包含未知参数的函数),一般情况下(包括伯努利,高斯),T(y)=y。a(
η
\eta
η)为一个和参数
η
\eta
η有关的函数。
η
\eta
η为一个标量,
η
T
\eta^T
ηT=
η
\eta
η,
η
T
\eta^T
ηT*T(y)=
η
T
\eta^T
ηT*y为两个标量相乘。在某些特殊(不常用~)分布中,
η
\eta
η可能为一个矩阵,T(y)是一个向量
由选定的指数分布族得到广义线性模型的过程需要三个假设:
1) y|x,
θ
\theta
θ~ExponentialFamily(
η
\eta
η),即对给定的输入x和参数
θ
\theta
θ时,y服从以
η
\eta
η为参数的指数分布族分布
2)对给定的输入x,输出预测函数
h
θ
(
x
)
=
E
(
T
(
y
)
∣
x
)
h_\theta(x)=E(T(y)|x)
hθ(x)=E(T(y)∣x),通常情况下T(y)=y,所以输出函数
h
θ
(
x
)
=
E
(
y
∣
x
)
h_\theta(x)=E(y|x)
hθ(x)=E(y∣x)即给定x时y的数学期望。
3)参数
η
\eta
η=
θ
T
x
\theta^Tx
θTx (若参数为向量不是标量,则
η
i
\eta_i
ηi=
θ
i
T
x
i
\theta_i^Tx_i
θiTxi)
由此,得到输出y的预测函数
h
θ
(
x
)
=
E
(
y
∣
x
)
h_\theta(x)=E(y|x)
hθ(x)=E(y∣x)
有两个较专业的术语:
g
(
η
)
=
E
(
y
∣
η
)
g(\eta)=E(y|\eta)
g(η)=E(y∣η)(即h_\theta(x)表达式中
θ
T
x
\theta^T x
θTx用
e
t
a
eta
eta来替换)称为正则响应函数,
g
−
1
(
η
)
g^{-1}(\eta)
g−1(η)称为正则关联函数。这是数学研究才需要了解的。
2、伯努利分布用广义线性模型来表示
伯努利分布
P
P
P{y=1;
θ
\theta
θ}=
θ
\theta
θ,可写成P{y;
θ
\theta
θ}=
θ
y
(
1
−
θ
)
1
−
y
\theta^y (1-\theta)^{1-y}
θy(1−θ)1−y=exp{
l
n
(
θ
y
(
1
−
θ
)
1
−
y
)
ln(\theta^y (1-\theta)^{1-y})
ln(θy(1−θ)1−y)}=exp{
y
∗
l
n
(
θ
)
+
(
1
−
y
)
∗
l
n
(
1
−
θ
)
y*ln(\theta)+(1-y)*ln(1-\theta)
y∗ln(θ)+(1−y)∗ln(1−θ)}=exp{
y
∗
l
n
(
θ
)
+
l
n
(
1
−
θ
)
−
y
∗
l
n
(
1
−
θ
)
y*ln(\theta)+ln(1-\theta)-y*ln(1-\theta)
y∗ln(θ)+ln(1−θ)−y∗ln(1−θ)}=exp{
y
∗
l
n
(
θ
1
−
θ
+
l
n
(
1
−
θ
)
)
y*ln(\frac{\theta}{1-\theta}+ln(1-\theta))
y∗ln(1−θθ+ln(1−θ))},y=0或1.此时,
l
n
(
θ
1
−
θ
)
ln(\frac{\theta}{1-\theta})
ln(1−θθ)为广义线性模型的参数
η
\eta
η。此时,求解参数
η
\eta
η可得
η
\eta
η=
l
n
(
θ
1
−
θ
)
ln(\frac{\theta}{1-\theta})
ln(1−θθ),即
θ
\theta
θ=
1
1
+
e
−
η
\frac{1}{1+e^{-\eta}}
1+e−η1。而
l
n
(
1
−
θ
)
ln(1-\theta)
ln(1−θ)为广义线性模型的函数
a
(
η
)
a(\eta)
a(η),可解得
a
(
η
)
=
−
l
n
(
1
−
θ
)
=
l
n
(
1
+
e
η
)
a(\eta)=-ln(1-\theta)=ln(1+e^\eta)
a(η)=−ln(1−θ)=ln(1+eη) T(y)=y,b(y)=1
由于伯努利分布的数学期望E(y)=
θ
\theta
θ,所以
h
θ
(
x
)
h_\theta(x)
hθ(x)=E(y|x)=
θ
\theta
θ=
1
1
+
e
−
η
\frac{1}{1+e^{-\eta}}
1+e−η1=
1
1
+
e
−
θ
T
x
\frac{1}{1+e^{-\theta^T x}}
1+e−θTx1这就是logistic函数的由来。伯努利分布的预测函数
h
θ
(
x
)
h_\theta(x)
hθ(x)=
1
1
+
e
−
θ
T
x
\frac{1}{1+e^{-\theta^T x}}
1+e−θTx1由此选定。
3、高斯分布用广义线性模型来表示
令高斯分布的方差
σ
2
\sigma^2
σ2=1
P
P
P{
y
;
μ
,
1
y;\mu,1
y;μ,1}=
1
2
π
e
−
(
y
−
μ
)
2
2
\frac{1}{\sqrt{2\pi}}e^{-\frac{(y-\mu)^2}{2}}
2π1e−2(y−μ)2 =
1
2
π
\frac{1}{\sqrt{2\pi}}
2π1 exp{
−
y
2
2
-\frac{y^2}{2 }
−2y2}exp{
μ
y
−
1
2
μ
2
\mu y-\frac{1}{2}\mu^2
μy−21μ2},T(y)=y,广义线性模型的参数
η
\eta
η=
η
T
\eta^T
ηT=
μ
\mu
μ。
a
(
η
)
=
1
2
μ
2
a(\eta)=\frac{1}{2}\mu^2
a(η)=21μ2,
b
(
y
)
=
1
2
π
b(y)=\frac{1}{\sqrt{2\pi}}
b(y)=2π1exp{
−
y
2
2
-\frac{y^2}{2}
−2y2}。
高斯分布的数学期望E(y|x)=
μ
\mu
μ,所以
h
θ
(
x
)
h_\theta(x)
hθ(x)=E(y|x)=
μ
\mu
μ=
η
\eta
η=
θ
T
x
\theta^T x
θTx。这就是线性回归预测函数
h
θ
(
x
)
h_\theta(x)
hθ(x)=
θ
T
x
\theta^T x
θTx的由来。
4、多项式分布(T(y)!=y)
多项式分布适用于有k个选择结果的分类问题。这类分布使用广泛,但T(y)!=y。
设k个选择结果的可能性分别为
ϕ
1
,
ϕ
2
.
.
.
ϕ
k
\phi_1,\phi_2...\phi_k
ϕ1,ϕ2...ϕk,由于
∑
i
=
1
i
=
k
ϕ
i
=
1
\sum_{i=1}^{i=k} \phi_i=1
∑i=1i=kϕi=1,所以只需要前k-1项。多项式分布的概率公式可以写为P{Y=i}=
ϕ
i
=
∏
i
=
1
i
=
k
−
1
ϕ
i
T
(
y
)
i
∗
ϕ
k
1
−
∑
j
=
1
j
=
k
−
1
T
(
y
)
j
\phi_i=\prod_{i=1}^{i=k-1} \phi_i^{T(y)_i}*\phi_k^{1-\sum_{j=1}^{j=k-1}T(y)_j}
ϕi=∏i=1i=k−1ϕiT(y)i∗ϕk1−∑j=1j=k−1T(y)j,其中,T(y)i=[0,0,0…1,0,0…0]T,是一个(k-1)维向量,第i项为1其余项为0,T(y)1=[1,0,0…0]T,T(y)k-1=[0,0,…0,0,1]T,T(y)k=[0,0,0…0,0]T定义为零向量。该概率公式连乘式中只有Y=i对应的
ϕ
i
\phi_i
ϕi上标为1,别的
ϕ
i
\phi_i
ϕi下标都为0。上式也可引入指示函数 I{True}=1,I{false}=0,则T(y)i=I(y=i),则P{Y=i}=
ϕ
i
=
∏
i
=
1
i
=
k
−
1
ϕ
i
I
(
y
=
k
)
∗
ϕ
k
1
−
∑
j
=
1
j
=
k
−
1
I
(
y
=
j
)
=
∏
i
=
1
i
=
k
ϕ
i
i
(
y
=
k
)
\phi_i=\prod_{i=1}^{i=k-1} \phi_i^{I(y=k)}*\phi_k^{1-\sum_{j=1}^{j=k-1}I(y=j)}=\prod_{i=1}^{i=k}\phi_i^{i(y=k)}
ϕi=∏i=1i=k−1ϕiI(y=k)∗ϕk1−∑j=1j=k−1I(y=j)=∏i=1i=kϕii(y=k)
把该分布写成指数分布族的形式,P{Y=i}=
∏
i
=
1
i
=
k
ϕ
i
I
(
y
=
k
)
\prod_{i=1}^{i=k} \phi_i^{I(y=k)}
∏i=1i=kϕiI(y=k)=exp{
T
(
y
)
1
∗
l
n
(
ϕ
1
)
+
T
(
y
)
2
∗
l
n
(
ϕ
2
)
+
.
.
.
T
(
y
)
k
−
1
∗
l
n
(
ϕ
k
−
1
)
+
(
1
−
∑
j
=
1
j
=
k
−
1
T
(
y
)
j
)
∗
l
n
(
ϕ
k
)
T(y)_1*ln(\phi_1)+T(y)_2*ln(\phi_2)+...T(y)_{k-1}*ln(\phi_{k-1})+(1-\sum_{j=1}^{j=k-1}T(y)_j)*ln(\phi_k)
T(y)1∗ln(ϕ1)+T(y)2∗ln(ϕ2)+...T(y)k−1∗ln(ϕk−1)+(1−∑j=1j=k−1T(y)j)∗ln(ϕk)}把
∑
j
=
1
j
=
k
−
1
T
(
y
)
j
∗
l
n
(
ϕ
i
)
\sum_{j=1}^{j=k-1}T(y)_j*ln(\phi_i)
∑j=1j=k−1T(y)j∗ln(ϕi)分别与前k-1项对应项相加,可以得到原式=exp{
T
(
y
)
1
∗
l
n
(
ϕ
1
ϕ
k
)
,
T
(
y
)
2
∗
l
n
(
ϕ
2
ϕ
k
)
,
.
.
.
T
(
y
)
k
−
1
∗
l
n
(
ϕ
k
−
1
ϕ
k
)
,
l
n
(
ϕ
k
)
T(y)_1*ln(\frac{\phi_1}{\phi_k}),T(y)_2*ln(\frac{\phi_2}{\phi_k}),...T(y)_{k-1}*ln(\frac{\phi_{k-1}}{\phi_k}),ln(\phi_k)
T(y)1∗ln(ϕkϕ1),T(y)2∗ln(ϕkϕ2),...T(y)k−1∗ln(ϕkϕk−1),ln(ϕk)}=exp{
η
T
T
(
y
)
y
+
l
n
(
ϕ
k
)
\eta^T T(y)_y +ln(\phi_k)
ηTT(y)y+ln(ϕk)}
由此得出该分布的T(y)=[0,0,0…1,0,0…0]^T(y=i对应的第i项为1其余为0,若y=k则为0向量),
η
=
[
l
n
(
ϕ
1
ϕ
k
)
,
l
n
(
ϕ
2
ϕ
k
)
,
.
.
.
,
l
n
(
ϕ
k
−
1
ϕ
k
)
]
T
\eta =[ln(\frac{\phi_1}{\phi_k}),ln(\frac{\phi_2}{\phi_k}),...,ln(\frac{\phi_{k-1}}{\phi_k})]^T
η=[ln(ϕkϕ1),ln(ϕkϕ2),...,ln(ϕkϕk−1)]T与T(y)一样为1个(k-1)维向量,,
a
(
η
)
=
−
l
n
(
ϕ
k
)
a(\eta)=-ln(\phi_k)
a(η)=−ln(ϕk),b(y)=1
η
i
=
l
n
(
ϕ
i
ϕ
k
)
(
i
属
于
[
1
,
k
−
1
]
)
,
e
η
i
=
ϕ
i
ϕ
k
\eta_i=ln(\frac{\phi_i}{\phi_k}) (i属于[1,k-1]),e^{\eta_i}=\frac{\phi_i}{\phi_k}
ηi=ln(ϕkϕi)(i属于[1,k−1]),eηi=ϕkϕi累加可得
∑
i
=
1
i
=
k
e
η
i
=
∑
i
=
1
i
=
k
ϕ
i
ϕ
k
=
1
ϕ
k
\sum_{i=1}^{i=k}e^{\eta_i}=\frac{\sum_{i=1}^{i=k}\phi_i}{\phi_k}=\frac{1}{\phi_k}
∑i=1i=keηi=ϕk∑i=1i=kϕi=ϕk1,由此可得
ϕ
k
=
1
∑
i
=
1
i
=
k
e
η
i
\phi_k=\frac{1}{\sum_{i=1}^{i=k}e^{\eta_i}}
ϕk=∑i=1i=keηi1。将解出的
ϕ
k
\phi_k
ϕk代入
e
η
i
=
ϕ
i
ϕ
k
e^{\eta_i}=\frac{\phi_i}{\phi_k}
eηi=ϕkϕi可得
ϕ
i
=
e
η
i
∑
i
=
1
i
=
k
e
η
i
\phi_i=\frac{e^{\eta_i}}{\sum_{i=1}^{i=k}e^{\eta_i}}
ϕi=∑i=1i=keηieηi。注意到多项式分布Y=i的数学期望就是
ϕ
i
\phi_i
ϕi,由此可以解得多项式分布的预测函数
h
θ
(
x
)
h_\theta(x)
hθ(x)=
E
(
T
(
y
)
∣
x
,
θ
)
E(T(y)|x,\theta)
E(T(y)∣x,θ)=[E(P{Y=1}),E(P{Y=2}),…E(P{Y=k-1})]T=
[
ϕ
1
,
ϕ
2
,
.
.
.
ϕ
k
−
1
]
T
=
[
e
η
1
∑
i
=
1
i
=
k
e
η
i
,
e
η
2
∑
i
=
1
i
=
k
e
η
i
.
.
.
e
η
k
−
1
∑
i
=
1
i
=
k
e
η
i
]
T
=
[
e
θ
1
T
x
∑
i
=
1
i
=
k
e
θ
i
T
x
,
e
θ
2
T
x
∑
i
=
1
i
=
k
e
θ
i
T
x
,
.
.
.
e
θ
k
−
1
T
x
∑
i
=
1
i
=
k
e
θ
i
T
x
]
T
[\phi_1,\phi_2,...\phi_{k-1}]^T=[\frac{e^{\eta_1}}{\sum_{i=1}^{i=k}e^{\eta_i}},\frac{e^{\eta_2}}{\sum_{i=1}^{i=k}e^{\eta_i}}...\frac{e^{\eta_{k-1}}}{\sum_{i=1}^{i=k}e^{\eta_i}}]^T=[\frac{e^{\theta_1^Tx}}{\sum_{i=1}^{i=k}e^{\theta_i^Tx}},\frac{e^{\theta_2^Tx}}{\sum_{i=1}^{i=k}e^{\theta_i^Tx}},...\frac{e^{\theta_{k-1}^Tx}}{\sum_{i=1}^{i=k}e^{\theta_i^Tx}}]^T
[ϕ1,ϕ2,...ϕk−1]T=[∑i=1i=keηieη1,∑i=1i=keηieη2...∑i=1i=keηieηk−1]T=[∑i=1i=keθiTxeθ1Tx,∑i=1i=keθiTxeθ2Tx,...∑i=1i=keθiTxeθk−1Tx]T
5、总结
当预测值服从高斯分布时,选用高斯分布的广义线性模型得出的 h θ ( x ) h_\theta(x) hθ(x)= θ T x \theta^T x θTx来进行预测,高维高斯分布即为参数 θ , η \theta,\eta θ,η与输入值x均为向量的情况;当预测值服从伯努利分布(0-1分布)时,选用伯努利分布的广义线性模型得出的 h θ ( x ) h_\theta(x) hθ(x)= 1 1 + e − θ T x \frac{1}{1+e^{-\theta^T x}} 1+e−θTx1进行预测。对于其他预测值分布形式已知或已建模的情形,(用于计数的泊松分布,用于估计时间发生时间的指数分布,用于计算时间间隔的Gamma分布等)均可以用广义线性模型推出合适的预测函数 h θ ( x ) h_\theta(x) hθ(x),以此编写代码。
2018.12.5