CS229 2017版作业2

课程视频地址:http://open.163.com/special/opencourse/machinelearning.html

课程主页:http://cs229.stanford.edu/

更具体的资料链接:https://www.jianshu.com/p/0a6ef31ff77a

作业地址:https://github.com/Doraemonzzz/CS229

参考资料:https://github.com/zyxue/stanford-cs229

这部分回顾2017版作业2。

1. Logistic Regression: Training stability

(a)在数据A上训练logistic regression model很快就收敛了,在数据B上训练logistic regression model无法收敛。

(b)观察后可以发现$\theta$的模长来越大,回顾logistic regression model

当$\theta$的模长很大时,$\theta^Tx$的模长很大,$g(\theta^T x) \to 0$,$g(z)’ = g(z)(1-g(z)) \to 0$,从而梯度会越来越小,训练会很慢。

之所以数据B发生这个现象而A没有发生这个现象,是因为数据A线性不可分,数据B线性可分。

由数据B线性可分可的

我们的目标函数为

要使得使得目标函数越小,只要$h_{\theta}(y^{(i)}x^{(i)})$越大即可,由于$y_i\theta^T x_i\ge 0$,所以$\theta$的模长越大,$y_i\theta^T x_i$就会越大,由梯度下降的性质可知,迭代之后会让$\theta​$的模长越来越大,就会发生上述现象。

而数据$A$不是线性可分的,所以存在$j​$,使得

所以算法不会让$\theta$的模长不停地增加。

(c)要解决上述问题,关键是不能让$\theta$的模长不停地增长,所以(iii),(v)是最好的方法。

(d)SVM不会发生这个问题,因为SVM是最大间隔分类器,即使可分,最大距离分类器也是唯一的,不会无限迭代下去。

而logistic回归实际上是在让函数间隔变大,所以会出现无法收敛的情形。

2.Model Calibration

(a)只要考虑两个分子即可,logistic回归的输出范围为$(0,1)$,题目中的$(a,b) = (0,1)$,所以

接下来证明这两项相等。

回顾损失函数

回顾课本的讲义可得

那么

其中

由$\theta​$的选择规则可知

这里有$n+1$个等式,注意$X^T$的第一行全为$1$,所以我们考虑第一个等式

由于$y^{(i)} \in \{0,1\}, h_θ(x^{(i)}) = P (y^{(i)} = 1|x^{(i)} ; θ)$,所以上式即为

从而

命题得证。

(b)考虑两个数据的数据集$x^{(1)},x^{(2)}$,不妨设$y^{(1)}=1,y^{(2)}=0$,如果

那么我们预测$y^{(1)} = 0,y^{(2)} = 1​$,准确率为$0​$,但是

所以perfectly calibrated无法推出perfect accuracy。

反之,如果

那么我们预测$y^{(1)} = 1,y^{(2)} = 0$,此时准确率为$1$,但是

所以perfect accuracy无法推出perfectly calibrated。

(c)设损失函数为

继续使用(a)的记号,那么

依旧考虑第一个等式

从而

3.Bayesian Logistic Regression and weight decay

回顾定义

由定义可知

因为

所以

从而

4.Constructing kernels

假设$K_i$对应的矩阵为$M_i$,$K$对应矩阵为$M$,由核函数的定义可知$M_i $为半正定阵。

(a)$K(x,z)= K_1(x,z) + K_2(x,z) $是核,因为

(b)$K(x,z)= K_1(x,z) -K_2(x,z) $不是核。取$K_2(x,z)=2K_1(x,z)$

(c)$K(x,z) = aK_1(x,z) ,a>0$是核

(d)$K(x,z) = -aK_1(x,z) ,a>0​$不是核

(e)$K(x,z) = K_1(x,z) K_2(x,z)$是核

因为$K_1,K_2$为核,所以设$K_1(x,z) = Φ_1(x) Φ^T_1(z),K_2=Φ_2(x) Φ^T_2(z)$。

记$\Phi(x)$是$Φ_1(x)Φ_2^T(x)$每一行拼接而成的向量,设$Φ_1(x),Φ_2(x)\in \mathbb R^n$,给出以下记号

那么

接着计算$\Phi(x) \Phi^T(x’) ​$,注意$\Phi^i(x)​$为行向量

所以$Φ(x)$对应的核为$K_1(x,x^{‘})K_2(x,x^{‘})$,从而$K(x,z) = K_1(x,z) K_2(x,z)$是核。

(f)$K(x,z) = f(x)f(z)$是核,因为符合定义。

(g)$K(x,z) = K_3(\phi(x), \phi(z))$是核,因为

(h)由(e)可知,如果$K_1$是核,那么$K_1^i (i\ge1, i\in N)$也是核,又由(a)(c)可得核函数的正系数的线性组合为核,所以$K(x,z)=p(K_1(x,z))$也是核。

5.Kernelizing the Perceptron

设这里的数据为$x^{(1)},…,x^{(m)}$

(a)根据更新公式

如果初始化$\theta^{(0)}=0​$,那么$\theta^{(i)}​$可以表示为$\phi(x^{(i)})​$的线性组合,从而

(b)计算$g({θ^{(i)}}^T \phi(x^{(i+1)}))​$

(c)由上述公式可知,第$i​$次我们只要更新$\beta_i ​$即可,更新公式如下

6.Spam classi fication

(a)(b)(c)

代码的注释比较详细,这里只说明一点,在朴素贝叶斯中我们需要计算:

题目中计算的对数概率,这是为了防止数值过小变为$0$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
# -*- coding: utf-8 -*-
"""
Created on Thu Mar 14 15:58:26 2019

@author: qinzhen
"""

import numpy as np
import matplotlib.pyplot as plt

def readMatrix(file):
fd = open(file, 'r')
#读取第一行
hdr = fd.readline()
#读取行列
rows, cols = [int(s) for s in fd.readline().strip().split()]
#读取单词
tokens = fd.readline().strip().split()
#构造空矩阵
matrix = np.zeros((rows, cols))
Y = []
#line为每行的元素
for i, line in enumerate(fd):
nums = [int(x) for x in line.strip().split()]
#第一个元素表示是否为垃圾邮件
Y.append(nums[0])
#将后续数据读入
kv = np.array(nums[1:])
#从第一个开始每两个累加
k = np.cumsum(kv[:-1:2])
#从第二个开始每隔一个取出
v = kv[1::2]
#这里应该是一种特殊的存储格式,我们直接使用即可
matrix[i, k] = v
return matrix, tokens, np.array(Y)

def nb_train(matrix, category):
#matrix(i,j)表示第i个邮件中第j个元素出现了几次
state = {}
#token的数量
N = matrix.shape[1]
#邮件数量
M = matrix.shape[0]
###################

#垃圾邮件的数量
y1 = matrix[category==1]
n1 = np.sum(y1)
#非垃圾邮件的数量
y0 = matrix[category==0]
n0 = np.sum(y0)

#P(y=1)
p1 = y1.shape[0] / M
#P(y=0)
p0 = y0.shape[0] / M
state[-1] = [p0, p1]

for i in range(N):
#找到第i个token
#垃圾邮件中第i个token出现的数量
s1 = matrix[category==1][:, i]
#拉普拉斯平滑
u1 = (s1.sum() + 1) / (n1 + N)
#非垃圾邮件中第i个token出现的数量
s0 = matrix[category==0][:, i]
#拉普拉斯平滑
u0 = (s0.sum() + 1) / (n0 + N)
#存入字典
state[i] = [u0, u1]

###################
return state

def nb_test(matrix, state):
output = np.zeros(matrix.shape[0])
###################
#邮件数量
M = matrix.shape[0]
#token的数量
N = matrix.shape[1]
for i in range(M):
#第i个邮件
vector = matrix[i]
s1 = np.log(state[-1][1])
s0 = np.log(state[-1][0])

for j in range(N):
#对第j个token的对数概率做累加
s1 += vector[j] * np.log(state[j][1])
s0 += vector[j] * np.log(state[j][0])
if s1 > s0:
output[i] = 1

###################
return output

def evaluate(output, label):
error = (output != label).sum() * 1. / len(output)
print('Error: %1.4f' % error)
return error

def nb(file):
trainMatrix, tokenlist, trainCategory = readMatrix(file)
testMatrix, tokenlist, testCategory = readMatrix('MATRIX.TEST')

state = nb_train(trainMatrix, trainCategory)
output = nb_test(testMatrix, state)

return evaluate(output, testCategory)

trainMatrix, tokenlist, trainCategory = readMatrix('MATRIX.TRAIN')
testMatrix, tokenlist, testCategory = readMatrix('MATRIX.TEST')

state = nb_train(trainMatrix, trainCategory)
output = nb_test(testMatrix, state)

evaluate(output, testCategory)

#problem b
b=[]
for i in range(1448):
b.append((i,np.log(state[i][1])-np.log(state[i][0])))

b.sort(key=lambda i:i[-1],reverse=True)
key = b[:5]

word = []
for i in key:
word.append(tokenlist[i[0]])

print(word)

#problem c
size = ['.50','.100','.200','.400','.800','.1400']
size1 = [50, 100, 200, 400, 800, 1400]
train = "MATRIX.TRAIN"
error = []
for i in size:
file = train+i
error.append(nb(file))

plt.plot(size, error)
plt.title("error VS szie")
1
2
3
4
5
6
7
Error: 0.0025
Error: 0.0238
Error: 0.0088
Error: 0.0088
Error: 0.0013
Error: 0.0000
Error: 0.0000

png

(d)这部分老师已经提供,但这里还是详细解释下。

首先这题需要计算高斯核矩阵,所以我们需要计算$[||x^{(i)}-x^{(j)}||^2]_{i,j}$,下面介绍向量化计算的方法:

假设

其中$x^{(i)} ,y^{(i)} \in \mathbb R^d$,现在的问题是如何高效计算矩阵$D \in \mathbb R^{m\times n}$,其中

首先对$D_{i,j}$进行处理

那么

特别的,这里$X=Y$,所以上述代码如下:

1
2
3
d1 = np.sum(matrix ** 2, axis=1).reshape(-1, 1)
d2 = d1.T
dist = d1 + d2 - 2 * matrix.dot(matrix.T)

带入高斯核的计算公式可得:

1
k = np.exp(- dist / (2 * (tau ** 2)))

老师给出的代码为:

1
2
3
gram = matrix.dot(matrix.T)
squared = np.sum(matrix*matrix, axis=1)
k = np.exp(-(squared.reshape((-1,1)) + squared.reshape((1,-1)) - 2 * gram) / (2 * (tau ** 2)))

为了解释剩余代码,首先介绍SVM对应的合页损失函数(参考统计学习方法113页):

其中

假设SVM对应的特征转换为$\phi (x)$,那么由课上的讨论可知,$w$的形式如下:

与之对应的损失函数为

由(1)可知,我们只要求出$\alpha= (\alpha_1,…,\alpha_N)^T$即可,所以可以关于$\alpha$做梯度随机梯度下降法,注意一个样本(不妨设为$x^{(i)}$)对应的损失函数为

记核矩阵为$K$,对上式关于$\alpha_k $求偏导可得

因此梯度为

所以计算梯度的过程如下,首先随机选个更新的样本(M为样本的数量):

1
i = int(np.random.rand() * M)

接着计算函数间隔(k为之前计算的核矩阵):

1
margin = category[i] * (k[i, :].dot(alpha))

然后计算正则项的梯度,注意这里我和老师给的计算式不大一样,老师给的式子如下:

1
grad = M * L * k[:, i] * alpha[i]

我的计算式如下:

1
grad = L * k.dot(alpha)

这里系数没有那么重要,实验结果表明两种方法的效果都尚可。

最后根据函数间隔是否大于$1$决定是否更新前一项的梯度:

1
2
if(margin < 1):
grad -= category[i] * k[:, i]

代码中的梯度除以了迭代次数$j​$的一个函数,这是为了让梯度逐渐变小:

1
alpha -= grad / ((np.sqrt(j+1)))

完整代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
# -*- coding: utf-8 -*-
"""
Created on Thu Mar 14 16:04:08 2019

@author: qinzhen
"""

import numpy as np
import matplotlib.pyplot as plt

tau = 8.

def readMatrix(file):
fd = open(file, 'r')
hdr = fd.readline()
rows, cols = [int(s) for s in fd.readline().strip().split()]
tokens = fd.readline().strip().split()
matrix = np.zeros((rows, cols))
Y = []
for i, line in enumerate(fd):
nums = [int(x) for x in line.strip().split()]
Y.append(nums[0])
kv = np.array(nums[1:])
k = np.cumsum(kv[:-1:2])
v = kv[1::2]
matrix[i, k] = v
#化为正负1
category = (np.array(Y) * 2) - 1
return matrix, tokens, category

def svm_train(matrix, category):
state = {}
M, N = matrix.shape
#####################
#大于0的化为1
matrix = 1.0 * (matrix > 0)
'''
#构造kernel矩阵
d1 = np.sum(matrix ** 2, axis=1).reshape(-1, 1)
d2 = d1.T
squared = matrix.dot(matrix.T)
dist = d1 + d2 - 2 * squared
k = np.exp(- dist / (2 * (tau ** 2)))
'''
gram = matrix.dot(matrix.T)
squared = np.sum(matrix*matrix, axis=1)
k = np.exp(-(squared.reshape((-1,1)) + squared.reshape((1,-1)) - 2 * gram) / (2 * (tau ** 2)))

#初始化向量
alpha = np.zeros(M)
#循环次数
n = 40
#系数
L = 1. / (64 * M)
#平均值
alpha_avg = np.zeros(M)

for j in range(n * M):
#随机取一个样本
i = int(np.random.rand() * M)
#计算函数间隔
margin = category[i] * (k[i, :].dot(alpha))
#grad = M * L * k[:, i] * alpha[i]
grad = L * k.dot(alpha)
if(margin < 1):
grad -= category[i] * k[:, i]
alpha -= grad / ((np.sqrt(j+1)))
alpha_avg += alpha

alpha_avg /= (n * M)

state['alpha'] = alpha
state['alpha_avg'] = alpha_avg
state['Xtrain'] = matrix
state['Sqtrain'] = squared

####################
return state


def svm_test(matrix, state):
M, N = matrix.shape
output = np.zeros(M)
###################
Xtrain = state['Xtrain']
Sqtrain = state['Sqtrain']
#大于0的化为1
matrix = 1.0 * (matrix > 0)
#做测试集的kernel
gram = matrix.dot(Xtrain.T)
squared = np.sum(matrix * matrix, axis=1)
k = np.exp(-(squared.reshape((-1, 1)) + Sqtrain.reshape((1, -1)) - 2 * gram) / (2 * (tau**2)))
#读取alpha
alpha_avg = state['alpha_avg']
#预测
pred = k.dot(alpha_avg)
output = np.sign(pred)

###################
return output

def evaluate(output, label):
error = (output != label).sum() * 1. / len(output)
print('Error: %1.4f' % error)
return error

def svm(file):
trainMatrix, tokenlist, trainCategory = readMatrix(file)
testMatrix, tokenlist, testCategory = readMatrix('MATRIX.TEST')

state = svm_train(trainMatrix, trainCategory)
output = svm_test(testMatrix, state)

return evaluate(output, testCategory)


trainMatrix, tokenlist, trainCategory = readMatrix('MATRIX.TRAIN.400')
testMatrix, tokenlist, testCategory = readMatrix('MATRIX.TEST')

state = svm_train(trainMatrix, trainCategory)
output = svm_test(testMatrix, state)

evaluate(output, testCategory)

size = ['.50','.100','.200','.400','.800','.1400']
size1 = [50, 100, 200, 400, 800, 1400]
train = "MATRIX.TRAIN"
error = []
for i in size:
file = train+i
error.append(svm(file))

plt.plot(size, error)
plt.title("error VS szie")
1
2
3
4
5
6
7
Error: 0.0025
Error: 0.0163
Error: 0.0075
Error: 0.0025
Error: 0.0025
Error: 0.0000
Error: 0.0000

png

(e)对比两图可以发现,和NB算法相比,SVM算法的效果更好。

本文标题:CS229 2017版作业2

文章作者:Doraemonzzz

发布时间:2019年03月14日 - 18:34:00

最后更新:2019年03月14日 - 18:39:11

原始链接:http://doraemonzzz.com/2019/03/14/CS229 2017版作业2/

许可协议: 署名-非商业性使用-禁止演绎 4.0 国际 转载请保留原文链接及作者。