机器学习基石作业03：二阶泰勒展开，特征转换，逻辑回归

本文总结机器学习基石的第三次作业，主要包括误差函数、特征转换、二阶泰勒展开，以及线性回归，逻辑回归的实现。

---

问题1和问题2：关于线性回归问题中$E_{in}$和$E_{out}$的理解。

关于含有噪声的目标$y=w_f^{T}x+ϵ$，其中的噪声$ϵ$均值为0，方差为${\sigma }^2$，且相互独立。根据PPT上的讲解可知，闭式解$w_{lin}$的$E_{in}$为：
$$E_{in}(w_{lin})=\frac{1}{N}||y-\hat{y}|{|}^2=\frac{1}{N}||(I-X{X}^{†})y|{|}^2=\frac{1}{N}||(I-H)y|{|}^2$$

从上图可知，$(I-H)y=(I-H)noise$（这是基于只有$f(x)$含有噪声，$x$不含噪声的前提），从而问题转换为$E_{in}(w_{lin})=\frac{1}{N}||(I-H)noise|{|}^2$。为了简化起见，令$A=(I-H)，B=noise\to ||AB|{|}^2=scaler$，从而可以获得下面的式子（其中用到$B{B}^T=scaler$）：
$$\begin{gathered} ||AB|{|}^2=trace((AB{)}^{T}AB)=trace(B{B}^T{A}^{T}A) \\ =B{B}^{T}trace(A^{T}A)=||noise|{|}^{2}trace(A^{T}A) \end{gathered}$$
根据$H$的性质（具体证明等见Q2）可得：
$$trace((I-H{)}^T(I-H))=trace(I-H)=N-(d+1)$$
所以，综上所述可得：
$$E_{in}(w_{lin})=(1-\frac{d+1}{N})||noise|{|}^2$$

---

当$\sigma =0.1,d=8$时，使得$E_{in}$的数学期望$\ge 0.008$的样本数是多少？（从选项中选择满足条件情况下最小的）

$E_{in}$的数学期望为：
$${\mathbb{E}}_D[E_{in}(w_{lin})]=(1-\frac{d+1}{N}){\sigma }^2$$
从而相当于$(1-9/N)\ast 0.01\ge 0.0008\to N=45$，从而选择N=100。

针对hat matrix$H=X(X^{T}X{)}^{-1}{X}^T$性质的探究，以下哪些性质是$H$所具有的？

结论：$H$是①对称 ②幂等性 ③半正定 ④有d+1个特征值为1。

证明：

• ①对称性：$H=(X(X^{T}X{)}^{-1}{X}^T{)}^T=X(X^{T}X{)}^{-1}{X}^T=H$（其中有用到$((AB{)}^{-1}{)}^T=((AB{)}^T{)}^{-1}$，这条式子有可逆作为先决条件）；
• ②幂等性：$H^2=X(X^{T}X{)}^{-1}{X}^{T}X(X^{T}X{)}^{-1}{X}^T=X(X^{T}X{)}^{-1}{X}^T=H$；
• ③半正定：假设存在特征值和特征向量$Hw=\lambda w$，则$\lambda w=Hw=H^{2}w=\lambda Hw={\lambda }^{2}w$，从而${\lambda }^{2}w=\lambda w\to \lambda =0or1$，所以全部特征值均$\ge 0$；
• ④存在d+1个特征值为1：$trace(H)=trace(X(X^{T}X{)}^{-1}{X}^T)=trace((X^{T}X{)}^{-1}{X}^{T}X)=trace(I_{d+1\times d+1})=d+1$，又根据$trace(H)=\sum {\lambda }_i$（该项需要半正定和对称性作为条件，具体证明可见AI圣经PRML），从而可知$\lambda =1$对应有$d+1$个。

这些结论可以用来证明$trace(I-H)=N-(d+1)$。

---

问题3-5：主要考察损失判据和随机梯度下降。


上述损失判据中，哪个是0/1判据$sign(w^{T}x)\ne y$的上界（ y ∈ { − 1 , + 1 } y\in \{-1,+1\} y∈{−1,+1}）？

通过上述图像易知，err3（本题D项，绘图存疑）即为0/1判据的上界。

---

上述损失判据中，哪个并不是处处可微分的？

由上图得，err2（本题D项）即在0处不可微分。

---

对上述损失判据计算SGD（忽略不可微分的情况），哪个损失判据恰好是PLA中采用的（即其求梯度函数恰好为PLA中更新参数时用到的）？

由第二节课件可知，PLA的参数更新方程为$w\leftarrow w+yxifsign(w^{T}x)\ne y$，可等价为$y{w}^{T}x>0,\nabla E=0,y{w}^{T}x<0,\nabla E=-yx$，从而推出$err(w)=max(0,-y{w}^{T}x)$

---

问题6-10：主要考查二元情况下的导数和二阶泰勒展开。

一阶导数：$\nabla f(x,y)=[\frac{\partial f}{\partial x},\frac{\partial f}{\partial y}{]}^T$
二阶导数：${\nabla }^{2}f(x,y)=[\frac{{\partial }^{2}f}{\partial x^2},\frac{{\partial }^{2}f}{\partial y\partial x};\frac{{\partial }^{2}f}{\partial y\partial x},\frac{{\partial }^{2}f}{\partial y^2}]$
二阶泰勒展开：
$$\begin{gathered} f(x+\Delta x,y+\Delta y)=f(x,y)+\Delta x\frac{\partial f(x,y)}{\partial x}+\Delta y\frac{\partial f(x,y)}{\partial y} \\ +\frac{1}{2!}[(\Delta x{)}^2\frac{{\partial }^{2}f(x,y)}{\partial x^2}+2\Delta x\Delta y\frac{{\partial }^{2}f(x,y)}{\partial x\partial y}+(\Delta y{)}^2\frac{{\partial }^{2}f(x,y)}{\partial y^2}] \end{gathered}$$
题设表达式：
$$E(u,v)=e^u+e^{2v}+e^{uv}+u^2-2uv+2v^2-3u-2v$$

$\nabla E(u,v)$在$(u,v)=(0,0)$处的值时多少？

根据一阶导数可得:
$$\begin{gathered} \frac{\partial E(u,v)}{\partial u}=e^u+v{e}^{uv}+2u-2v-3 \\ \frac{\partial E(u,v)}{\partial v}=2e^{2v}+u{e}^{uv}-2u+4v-2 \end{gathered}$$
将$(u,v)=(0,0)$代入可得 $\nabla E(0,0)=(-2,0)$。

---

根据梯度下降算法（如下式所示），对参数进行迭代更新，求$\eta =0.01,(u_0,v_0)=(0,0)$经过五次迭代后的结果$(u_5,v_5)$和$E(u_5,v_5)$？
$$(u_{t+1},v_{t+1})=(u_t,v_t)-\eta \nabla E(u_t,v_t)$$
由程序推出：$u=0.0941,v=0.0018,E=2.825$。

---

如果采用二阶泰勒展开${\hat{E}}_2(\Delta u,\Delta v)$来近似$E(u+\Delta u,v+\Delta v)$，求下述表达式中的参数在$(u,v)=(0,0)$处的值？
$${\hat{E}}_2(\Delta u,\Delta v)=b_{uu}(\Delta u{)}^2+b_{vv}(\Delta v{)}^2+b_{uv}(\Delta u)(\Delta v)+b_u\Delta u+b_v\Delta v+b$$
根据二阶导数的情况：
$$\begin{gathered} \frac{{\partial }^{2}E}{\partial u^2}=e^u+v^2{e}^{uv}+2 \\ \frac{{\partial }^{2}E}{\partial v^2}=4e^{2v}+u^2{e}^{uv}+4 \\ \frac{{\partial }^{2}E}{\partial u\partial v}=e^{uv}+vu{e}^{uv}-2 \end{gathered}$$
将上式代入二阶泰勒展开，推出：$(1.5,4,-1,-2,0,3)$。

---

将黑塞/海森矩阵（Hessian Matrix）表示为${\nabla }^{2}E(u,v)$，并假设该矩阵是正定的。以下哪个时最佳的$(\Delta u,\Delta v)$使得${\hat{E}}_2$取到最小值？（这个方向称之为Newton Direction）

找使得${\hat{E}}_2$最小的$(\Delta u,\Delta v)$，可以直接通过求导获得。分别对$\Delta u$和$\Delta v$求导，结果如下：
$$\begin{gathered} \frac{\partial E_2}{\partial \Delta u}=\frac{\partial E}{\partial u}+\Delta u\frac{{\partial }^{2}E}{\partial u^2}+\Delta v\frac{{\partial }^{2}E}{\partial u\partial v}=0 \\ \frac{\partial E_2}{\partial \Delta v}=\frac{\partial E}{\partial v}+\Delta v\frac{{\partial }^{2}E}{\partial v^2}+\Delta u\frac{{\partial }^{2}E}{\partial u\partial v}=0 \end{gathered}$$
联立上述两式，并将$(\Delta u,\Delta v)$以向量形式提取出来可以化简为：
$$(\Delta u,\Delta v{)}^T=-({\nabla }^{2}E{)}^{-1}\nabla E$$

---

对$(u_0,v_0)$利用Newton Direction（无$\eta$）进行参数更新，五轮更新后$(u_5,v_5)$和$E$ 结果为？

通过简单的程序可得：$u=0.6118,v=0.0705,E=2.3608$。

从该结果可以看出，Newton Direction更新速度更快，但是代价是求Hessian矩阵引入的复杂性。

---

问题11-12：关于特征转换的问题。

考虑二维空间上的6个点，采用含有二次函数，线性函数的hypotheses集合，最多能shatter其中的几个点。shatter的意思是说将输入样本点完全二分。看下图：

使用上述六条线可以将六个样本点shatter。

---

假设转换之前预先“偷看”了所有$N$个数据，并定义一种特殊的特征转换，将$x\in {\mathbb{R}}^d\to z\in {\mathbb{R}}^N$
$$(\Phi (x){)}_n=z_n=[x=x_n]$$
这题的关键在于理解这种“奇葩”的特征转换，举个例子，如第1个数据$x_1$，根据上述规则则变为$[1,0,...,0{]}^T$，（矩阵大小$N\times 1$）就是将第几个数对应的行置为1，其他行均为0。显然，不管多少数，其转换后的向量是两两正交的，正交的向量一定是线性无关的，因此均可以被shatter，所以$d_{vc}(H_{\Phi })=\infty$。所以答案为C。

---

问题13-15：主要关于线性回归问题和特征转换。

数据产生：数据集大小$N=1000$，且$\mathcal{X}=[-1,1]\times [-1,1]$，每个数据的$x$均等概率的从$\mathcal{X}$中提取。而对应的$y$则根据$f(x_1,x_2)=sign(x_1^2+x_2^2-0.6)$来确定，且对数据集中的$10\%$的数据的$y$进行反转（相当于添加噪声）。

先对线性回归算法进行简单的说明：

函数集：$y=w^{T}x$

损失函数：$E_{in}(w)=\frac{1}{N}{\sum }_{n=1}^N(w^T{x}_n-y_n{)}^2$

梯度：$\nabla E_{in}(w)=\frac{2}{N}(X^{T}Xw-X^{T}y)$

“目的”：寻找$w$使得损失函数最小

Linear Regression

①获得数据$(x_1,y_1),...,(x_N,y_N)$
②采用闭式解公式求出最佳$w$：$w_{lin}=(X^{T}X{)}^{-1}{X}^{T}y$
③返回$w_{lin}$

如果还有预测过程，直接$\hat{y}=w_{lin}^{T}x$

---

不进行特征转换，只采用特征$(1,x_1,x_2)$，利用Linear Regression获得最佳的$w_{lin}$。将其直接运用到分类问题上面（利用$sign(w^{T}x)$），在利用$0/1$判据来衡量训练样本误差$E_{in}$。进行1000次实验，取误差的平均。

Ein: 0.5036

通过上面结果可知，直接利用Linear Regression(利用square error)再运用到分类问题上结果很差！

---

问题14-15：将数据的特征进行转换，转换为$(1,x_1,x_2,x_1{x}_2,x_1^2,x_2^2)$这6项，再利用Linear Regression获得最佳的$w_{lin}$，求该$w_{lin}$以及将其运用到测试集上的测试误差$E_{out}$（衡量方式与Q13相同）。

theta:[[-1.01626639 0.07325707 0.02834912 -0.0155599 1.63387468 1.52477431]]
选项中最接近的为：
$$g(x_1,x_2)=sign(-1-0.05x_1+0.08x_2+0.13x_1{x}_2+1.5x_1^2+1.5x_2^2)$$

---

在14题得到的最优w的基础上，产生1000个测试样本，计算误差。重复1000次求平均。

Eout: 0.125225

---

问题16-17：关于多类别logistics regression问题。针对K类别分类问题，我们定义输出空间 Y = { 1 , 2 , . . . , K } \mathcal{Y}=\{1,2,...,K\} Y={1,2,...,K}，MLR的函数集可以视为由一系列(K个)权值向量$(w_1,...,w_K)$构成，其中每个权值向量均为$d+1$维。每种假设函数可以表示为：
$$h_y(x)=\frac{exp(w_y^{T}x)}{{\sum }_{i=1}^{K}exp(w_i^{T}x)}$$
且可以用来近似潜在的目标分布函数$P(y|x)$。MLR的“目标”就是从假设函数集中寻找使得似然函数最大的额假设函数。

类似Lecture10中最小化$-log(likelihood)$一样，推导$E_{in}(w_1,...,w_K)$。

采用同样的处理方式
$$max\frac{1}{N}\prod _{i=1}^N{h}_y(x)\to min-\frac{1}{N}\sum _{i=1}^{N}log(h_y(x))$$
将MLR的假设函数代入上式并化简可得：
$$\frac{1}{N}\sum _{n=1}^N(ln(\sum _{i=1}^{K}exp(w_i^T{x}_n))-w_{y_n}^T{x}_n)$$

---

针对上述的$E_{in}$，它的一阶导数$\nabla E_{in}$可以表示为$(\frac{\partial E_{in}}{\partial w_1},\frac{\partial E_{in}}{\partial w_2,},...,\frac{\partial E_{in}}{\partial w_K})$，求$\frac{\partial E_{in}}{\partial w_i}$。

直接对A16的答案的式子进行求导，就可以得到下式：
$$\frac{1}{N}\sum _{n=1}^N((h_i(x_n)-[y_n=i]x_n)$$

---

问题18-20：关于logistic regression实现的问题。

首先看一下算法：

函数集：$s={\sum }_{i=0}^d{w}_i{x}_i$，$h(x)=\theta (s)=\frac{1}{1+e^{-s}}$

损失函数：$E_{in}(w)=\frac{1}{N}{\sum }_{i=1}^{N}ln(1+exp(-y_n{w}^T{x}_n))$

梯度：$\nabla E_{in}=\frac{1}{N}{\sum }_{i=1}^N\theta (-y_n{w}^T{x}_n)(-y_n{x}_n)$

目标：寻找一个最佳假设函数使得损失函数最小

（注：$h(x)$来近似$P(y|x)$上述的损失函数通过cross-entropy可推导出来）

Logistic Regression：

初始化$w$
For t=0,1,…
① 计算$\nabla E_{in}(w)$
② 更新参数：$w\leftarrow w-\eta \nabla E_{in}(w)$
返回$w$

（上述$\eta$可以视为一个超参数，可以通过cross-validation来确定）

---

针对$\eta =0.001,T=2000$的情况，采用梯度下降法获得$w$后，在测试集上的错误率是多少？（利用0/1判据）

Ein = 0.466；Eout = 0.475。

---

针对$\eta =0.01,T=2000$的情况，采用梯度下降法获得$w$后，在测试集上的错误率是多少？（利用0/1判据）

Ein = 0.197；Eout = 0.22。

---

针对$\eta =0.001,T=2000$的情况，采用随机梯度下降法(此处采用按顺序每次选择元素，更通常的做法是随机选择元素)获得$w$后，在测试集上的错误率是多少？（利用0/1判据）

Ein = 0.464；Eout = 0.473。

---

import numpy as np
import pandas as pd
import math
import scipy.linalg as lin # 该模块包含线性代数的函数
import matplotlib.pyplot as plt

# Q3：4种不同的误差衡量和0/1误差
plt.figure(dpi=150)
plt.style.use('science')
z = np.arange(-2, 2, 0.01)

err0 = z.copy()
err0[err0 >= 0] = 0
err0[err0 < 0] = 1

err1 = 1-z.copy()
err1[err1 < 0] = 0

err2 = -z.copy()
err2[err2 < 0] = 0

err3 = 1-z.copy()
err3[err3 < 0] = 0
err3 = np.power(err3, 2)

err4 = 1/2 * (np.exp(-z))

plt.plot(z, err0, label='err0/1')
plt.plot(z, err1, label='err1')
plt.plot(z, err2, label='err2')
plt.plot(z, err3, label='err3')
plt.plot(z, err4, label='err4')

plt.legend()
plt.show()

# Q7
u = 0
v = 0
eta = 0.01
for i in range(5):
    du = math.exp(u) + v*math.exp(u*v) + 2*u-2*v-3
    dv = 2*math.exp(2*v) + u*math.exp(u*v) - 2*u + 4*v-2
    u -= eta*du
    v -= eta*dv
print('u=', np.round(u,decimals=4))
print('v=', np.round(v,decimals=4))
E = math.exp(u)+math.exp(2*v)+math.exp(u*v)+u**2-2*u*v+2*v**2-3*u-2*v
print(np.round(E,decimals=4))

# Q10
u = 0; v = 0
uv = np.array([[0], [0]])

for i in range(5):
    du = math.exp(u) + v*math.exp(u*v) + 2*u - 2*v - 3
    dv = 2*math.exp(2*v) + u*math.exp(u*v) - 2*u + 4*v-2

    du2 = math.exp(u) + v**2*math.exp(u*v) + 2
    dv2 = 4*math.exp(2*v) + u**2*math.exp(u*v) + 4
    dudv = math.exp(u*v) + v*u*math.exp(u*v) - 2
    
    ddE = np.array([[du2, dudv], [dudv, dv2]])
    dE = np.array([[du], [dv]])
    uv = uv-lin.inv(ddE).dot(dE) # linalg.inv()：矩阵求逆
    u = uv[0, 0]
    v = uv[1, 0]
print('u and v: ',uv.T)
E = math.exp(u) + math.exp(2*v) + math.exp(u*v) + u**2 - 2*u*v + 2*v**2 - 3*u -2*v
print('E: ',E)

# 数据生成函数
def generate_data(num):
    axeX = np.random.uniform(-1, 1, num) # 在[-1,1)内随机采样
    axeY = np.random.uniform(-1, 1, num)
    
    # np.c_：按【列】连接两个矩阵，要求行数相等。np.r_：按【行】连接两个矩阵，要求列数相等；
    Xtemp = np.c_[axeX, axeY]
    X = np.c_[np.ones((num, 1)), Xtemp]
    
    Ytemp = np.sign(np.power(axeX, 2)+np.power(axeY, 2)-0.6)
    Ytemp[Ytemp == 0] = -1
    
    pos = np.random.permutation(num)
    Ytemp[pos[0: round(0.1*num)]] *= -1
    
    Y = Ytemp.reshape((num, 1))
    
    return X, Y

totalerr = 0
for i in range(1000):
    X, Y = generate_data(1000)
    
    theta = lin.pinv(X.T.dot(X)).dot(X.T).dot(Y) # linalg.pinv矩阵伪逆
    
    ypred = np.sign(X.dot(theta))
    
    err = np.sum(ypred!=Y)/1000
    
    totalerr += err
    
print('Ein: ', totalerr/1000)

# 特征转换函数
def transform(X):
    row, col = X.shape
    Xback = np.zeros((row, 6))
    Xback[:, 0:col] = X
    Xback[:, col] = X[:, 1]*X[:, 2]
    Xback[:, col+1] = X[:, 1]**2
    Xback[:, col+2] = X[:, 2]**2
    return Xback

# Q14
totalerr = 0
for i in range(1000):
    X, Y = generate_data(1000)
    Xtran = transform(X)
    
    theta = lin.pinv(Xtran.T.dot(Xtran)).dot(Xtran.T).dot(Y)
    
    Xtest, Ytest = generate_data(1000)
    
    Xback = transform(Xtest)
    
    ypred = np.sign(Xback.dot(theta))
    
    err = np.sum(ypred!=Ytest)/1000
    
    totalerr += err
    
print('theta: ', theta.T)
print('Ein: ', totalerr/1000)

# sigmoid函数
def sigmoid(z):
    zback = 1/(1+np.exp(-1*z))
    return zback


# Logistic Regression
def logistic_regression(X, Y, eta, numiter, flag=0):
    row, col = X.shape
    theta = np.zeros((col, 1))
    num = 0
    for i in range(numiter):
        if flag == 0:
            derr = (-1*X*Y).T.dot(sigmoid(-1*X.dot(theta)*Y))/row
        else:
            if num >= row:
                num = 0
            derr = -Y[num, 0]*X[num: num+1, :].T*sigmoid(-1*X[num, :].dot(theta)[0]*Y[num, 0])
            num += 1
        theta -= eta*derr
    return theta

# 导入数据函数
def load_data(filename):
    data = pd.read_csv(filename, sep='\s+', header=None)
    col, row = data.shape
    X = np.c_[np.ones((col, 1)), data.iloc[:, 0: row-1]]
    Y = data.iloc[:, row-1:row].values
    return X, Y

# 误差计算函数
def mistake(X, Y, theta):
    yhat = X.dot(theta)
    yhat[yhat > 0] = 1
    yhat[yhat <= 0] = -1
    err = np.sum(yhat != Y)/len(Y)
    return err

X, Y = load_data('hw3_train.dat')
testX, testy = load_data('hw3_test.dat')

# Q18
eta = 0.001; T = 2000; flag = 0
theta = logistic_regression(X, Y, eta, T, flag)
errin = mistake(X, Y, theta)
errout = mistake(testX, testy, theta)
print(f'Ein = {errin};', 'Eout = ', errout)

# Q19
eta = 0.01; T = 2000; flag = 0
theta = logistic_regression(X, Y, eta, T, flag)
errin = mistake(X, Y, theta)
errout = mistake(testX, testy, theta)
print(f'Ein = {errin};', 'Eout = ', errout)

# Q20
eta = 0.001; T = 2000; flag = 1
theta = logistic_regression(X, Y, eta, T, flag)
errin = mistake(X, Y, theta)
errout = mistake(testX, testy, theta)
print(f'Ein = {errin};', 'Eout = ', errout)

---

参考：
机器学习基石作业3：https://github.com/AceCoooool/MLF-MLT