CS229 老版作业4

课程视频地址: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

这部分回顾老版作业4。

1. EM for supervised learning

(a)由定义可得

概率似然函数为

对数似然函数为

关于$\theta_i $求梯度可得

令上式为$0$可得

关于$\phi $求梯度可得

接着求Hessian矩阵,注意上述梯度第$s$个分量为$\sum_{j=1}^m\big(z^{(j)}-g(\phi^T x^{(j)})\big) x^{(j)}_s $,求二阶偏导数:

那么关于$\phi $的梯度为

Hessian矩阵为

(b)

概率似然函数为

对数似然函数为

为了使用EM算法,我们首先求后验概率

事实上,我们有

上述每一项都可以根据定义计算得到。

为了方便叙述,记

那么由EM算法,接下来我们需要最大化

记上式为$l_1$,关于$\theta_i $求梯度可得

其中最后一步是因为$z^{(j)}=k$。令上式为$0$可得

关于$\phi $求梯度可得

最后一步利用到了

那么关于$\phi $的梯度以及Hessian矩阵和(a)的形式一样:

关于$\phi $的梯度为

Hessian矩阵为

2. Factor Analysis and PCA

(a)我们可以等价的定义

由$z \sim \mathcal N(0,I)$我们知道$\mathbb E[z]=0$。所以,我们有

将上述结果合并,我们有

接着求协方差。因为$z \sim \mathcal N(0,I)$,我们可以轻松得到$\Sigma_{zz}=\text{Cov}(z)=I$。并且,我们有

类似的,我们可以按如下方式计算出$\Sigma_{xx}$

回顾之前的结论

所以

(b)不难发现$x$的边际分布为

所以对数似然函数为

接着使用EM算法,E步骤很简单,只要取

由之前讨论,不难发现

接着,我们需要最大化

弃不依赖参数$U$的项($Q_i$虽然包含$U$,但是在执行更新时是固定的),可以发现我们需要最大化:

让我们关于$U$最大化上式,这里需要利用

注意到只有最后一项依赖$U$,求梯度可得

令上式为$0$可得

其中

(c)如果$\sigma^2 \to 0$,那么

如果记

所以E步骤可以化为

M步骤可以化为

如果在EM步骤后$U$不变,那么我们有

注意这里$\frac {U^TU} {U^TX^TXU}$为标量,记其为$\frac 1 \lambda $,那么

所以结论成立。

3. PCA and ICA for Natural Images

这里要对讲义中ICA的更新公式稍作变形,首先回顾之前的公式:

对于训练样本$x^{(i)}$,随机梯度下降法为

那么

所以可以利用该矩阵得到

对应代码为

1
2
WX = W.dot(X1.T)
G = 1 - 2 * sigmoid(WX)

现在考虑整体梯度

化简可得上式为

对应代码为:

1
grad = G + batch * inv(W.T)

这里我将matlab代码改写为python,下面是代码:

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
import numpy as np
from numpy.linalg import eigh, inv, svd, norm
import matplotlib.pyplot as plt
from matplotlib.image import imread

patch_size = 16

X_ica = np.zeros((patch_size*patch_size, 40000))
idx = 0
for s in range(1, 5):
#转换为浮点型
image = imread("./images/{}.jpg".format(s)).astype(np.float32)
#获得数据维度
a, b, c = np.shape(image)
y = a
x = b * c
#注意这里order要选择F
image = image.reshape(y, x, order="F")
for i in range(1, y//patch_size+1):
for j in range(1, x//patch_size+1):
patch = image[(i-1)*patch_size: i*patch_size, (j-1)*patch_size: j*patch_size]
X_ica[:, idx] = patch.reshape(-1, 1, order="F").flatten()
idx += 1

X_ica = X_ica[:, 0: idx]

#正定矩阵分解
W = 1 / X_ica.shape[1] * X_ica.dot(X_ica.T)
w, v = eigh(W)
w = np.diag(1 / np.sqrt(w))
W_z = v.dot(w).dot(v.T)

X_ica = X_ica - np.mean(X_ica, axis=1).reshape(-1, 1)
X_pca = X_ica

X_ica = 2 * W_z.dot(X_ica)
X_pca = X_pca / np.std(X_pca, axis=1).reshape(-1, 1)

之前为预处理函数,可以忽略细节。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#### PCA
def pca(X):
U, S, V = svd(X.dot(X.T))

return U

def plot_pca_filters(U):
n = (patch_size + 1) * patch_size - 1
big_filters = np.min(U) * np.ones((n, n))
for i in range(patch_size):
for j in range(patch_size):
big_filters[i*(patch_size+1): (i+1)*(patch_size+1)-1, j*(patch_size+1): (j+1)*(patch_size+1)-1] = U[:, i*patch_size+j].reshape(patch_size, patch_size)

plt.imshow(big_filters)
plt.gray()
plt.show()

U = pca(X_pca)
plot_pca_filters(U)

png

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
#### ICA
def sigmoid(x):
return 1 / (1 + np.exp(-x))

X = X_ica.copy()

X = X.T
n, d = X.shape
W = np.eye(d)
batch = 100
alpha = 0.0005
eps = 1e-1

for i in range(10):
print("第{}轮".format(i+1))
np.random.permutation(X)
for j in range(n // batch):
X1 = X[j * batch: (j+1) * batch, :]
WX = W.dot(X1.T)
grad = (1 - 2 * sigmoid(WX)).dot(X1) + batch * inv(W.T)
W += alpha * grad


def plot_ica_filters(W):
F = W.dot(W_z)
norms = norm(F, axis=1)
idxs = np.argsort(norms)
norms = np.sort(norms)
n = (patch_size + 1) * patch_size - 1
big_filters = np.min(W) * np.ones((n, n))
for i in range(patch_size):
for j in range(patch_size):
temp = W[idxs[i*patch_size+j], :].reshape(patch_size, patch_size)
big_filters[i*(patch_size+1): (i+1)*(patch_size+1)-1, j*(patch_size+1): (j+1)*(patch_size+1)-1] = temp

plt.imshow(big_filters)
plt.gray()
plt.show()

plot_ica_filters(W)
1
2
3
4
5
6
7
8
9
10
第1轮
第2轮
第3轮
第4轮
第5轮
第6轮
第7轮
第8轮
第9轮
第10轮

png

4. Convergence of Policy Iteration

(a)利用定义即可:

(b)注意由定义我们有

所以

因此

(c)由$\pi’$的定义可得

所以

由(a)可得

所以

对右边不断作用$B^{\pi’}$,然后结合(b),我们得到

(d)因为状态和行动有限,所以

的取值有限,由(c)可得策略迭代的算法会导致$V^\pi(s)$非降,即

所以最终必然会达到

所以

带入$V^{\pi}(s)$的定义可得

5. Reinforcement Learning: The Mountain Car

题目的公式有误,正确的如下:

这里没有将代码改写为python,matlab代码参考了标准答案,结果如下:

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
function [q, steps_per_episode] = qlearning(episodes)

% set up parameters and initialize q values
alpha = 0.05;
gamma = 0.99;
num_states = 100;
num_actions = 2;
actions = [-1, 1];
q = zeros(num_states, num_actions);
steps_per_episode = zeros(1, episodes);

for i=1:episodes,
[x, s, absorb] = mountain_car([0.0 -pi/6], 0);
%%% YOUR CODE HERE
%%% 找到第一步的动作对应的最大值和索引
[maxq, a] = max(q(s, :));
%%% 如果相同则随机
if q(s, 1) == q(s, 2)
a = ceil(rand * num_actions);
end
%%% 更新步数
steps = 0;

%%% 如果未吸收
while(~absorb)
%%% 找到下一步的位置
[x, sn, absorb] = mountain_car(x, actions(a));
%%% 奖励
reward = - double(~ absorb);
%%% 找到动作对应的最大值和索引
[maxq, an] = max(q(sn, :));
%%% 如果相同则随机
if q(s, 1) == q(s, 2)
an = ceil(rand * num_actions);
end
%%% 找到最大的行动
q(s, a) = (1 - alpha) * q(s, a) + alpha * (reward + gamma * maxq);
%%% 更新状态
a = an;
s = sn;
steps = steps + 1;
end
%%% 记录步数
steps_per_episode(i) = steps;
end

本文标题:CS229 老版作业4

文章作者:Doraemonzzz

发布时间:2019年03月29日 - 20:04:00

最后更新:2019年03月29日 - 20:05:29

原始链接:http://doraemonzzz.com/2019/03/29/CS229 老版作业4/

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