课程视频地址：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)}$，随机梯度下降法为

记

那么

所以可以利用该矩阵得到

对应代码为

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

现在考虑整体梯度

化简可得上式为

对应代码为：

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

这里我将matlab代码改写为python，下面是代码：

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)

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

#### 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

#### 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轮

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代码参考了标准答案，结果如下：

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