mac_test

发表于 2019-08-25

entropy

通过概率

发表于 2019-04-27 | 分类于 code

问题

一共有100道题，第i道题做对的概率为p[i]，至少做对60道题为及格，求及格的概率？

代码

def passProbability(p,thresh=60):
    """
    f[i][j]: 前i道题至少对j道的概率
    :param p: probability list [p[0],...,p[N-1]]
    :param thresh: 及格线做对题目个数
    :return:
    """
    N = len(p)
    f = [[0]*(N+1) for _ in range(N+1)]
    # 0道题对0道
    f[0][0] = 1
    # i道题对0道
    for i in range(1,N+1):
        f[i][0] = f[i-1][0]*(1-p[i-1])
    # i道题对j道
    for i in range(1,N+1):
        for j in range(1,i+1):    # 1 =<j <=i 
            f[i][j] = f[i-1][j]*(1-p[i-1])+f[i-1][j-1]*p[i-1]
    # 及格概率
    result = 0
    for j in range(thresh,N+1):
        result += f[N][j]
    return result

及格概率与每道题做对概率的关系(假设每道题做对概率一样)

import  matplotlib.pyplot as plt
xs = [x*0.01 for x in range(1,100)]
ys = [passProbability([x]*100) for x in xs]
plt.figure()
plt.plot(xs,ys)
plt.show()

结果

分析

转移方程

设$f[i][j]$是前i道题做对j道的概率，分如下两种情况讨论

当$j\ge1$时
- 前i道对j-1道，最后一道对
- 前i-1道对j道，最后一道错
当$j=0$时：前i-1道对0道，最后一道错

从而有

$f[i][j] = \begin{cases} f[i][j-1]*p[i]+f[i-1][j]*(1-p[i]) \quad j\ge 1 \\ f[i-1][0]*(1-p[i])\quad j=0 \end{cases}$

边界条件

$f[0][0]=1$

$f[0][j]=0 (j\ge1)$

及格概率

100道题对至少60道为及格，所以及格概率为$\sum_{j=1}^{60} f[100][j]$

卡特兰数-Catalan number

发表于 2019-04-12 | 分类于 code

问题

栈是一种先进后出（FILO，First In Last Out）的数据结构。当栈的进栈顺序为$1,2,3,…,n$ 时，问有多少种不同的出栈序列？

分析

设$C(n)$是$n$个数所有合法出栈序列的个数。

思路1：递归

以某个基准数在出栈序列中的顺序来讨论，基准数用$k$表示。$C(n)$可以分解为：

$k$在最后一个位置。那么$k$先进栈，最后出栈；其他$n-1$个元素$12,k-1,k+1…,n$(共$n-1$个元素)的出栈个数为$C(n-1)$
$k$在倒数第二个位置。那么有一个元素比$k$先进栈，共$C(1)$种可能；其他$n-2$个元素的出栈个数为$C(n-2)$。由组合数的乘法原理，得到总的出栈个数为$C(1)*C(n-2)$
$\cdots$
$k$在第一个位置。那么$k$先进栈，然后马上出栈；其他$n-1$个元素的出栈个数为$C(n-1)$。总的出栈个数为$C(n-1)$。

令$C(0)=1$，从有

$C(n) = \sum_{i=0}^{n-1} C(i)C(n-i)$

参考母函数推卡特兰数通项公式^note1，可得

$C(n) = \frac{1}{n+1} \begin{pmatrix} 2n \\ n \end{pmatrix}$

思路2：排列组合

用$+1$表示进栈，$-1$表示出栈。那么出栈序列$3,2,1$(假设序列为$1,2,3$)的进出栈顺序为

$+1,+1,+1,-1,-1,-1$

对于$n$个数的序列，共有$n$次进栈($+1$)和$n$次出栈$-1$.

对$+1,-1$这$2n$个数字，对其做排列，我们的目标是找到所有合法的排列(对应合法的出栈顺序)。

总的排列个数为$\begin{pmatrix}
2n \\
n
\end{pmatrix}$
不合法的个数

不合的排列一定会出现这种情况：从左往右若干项和为$-1$, 该项一定在奇数位置。

共$n$个$+1$，$n$个$-1$

构造映射：将前$2i+1$项取反(第$2i+1$项为首次出现累加和为$-1$的位置), 结果为

共$n+1$个$+1$, $n-1$个$-1$

该映射为一一映射，所以不合法个数与共$n+1$个$+1$, $n-1$个$-1$的全排列个数一样，即，$\begin{pmatrix}
2n \\
n+1
\end{pmatrix}$或者$\begin{pmatrix}
2n \\
n-1
\end{pmatrix}$.

经过简单的计算，有

$C(n)= \begin{pmatrix} 2n \\ n \end{pmatrix}- \begin{pmatrix} 2n \\ n+1 \end{pmatrix}= \frac{1}{n+1} \begin{pmatrix} 2n \\ n \end{pmatrix}$

一些类似的问题

合法括号表达式个数

计算$n$对括号形成的合法括号表达式的个数。

合法：$(())(), ()()$

不合法：$((), ())))$

分析：将$($看做进栈，将$)$看做出栈。

答案： $C(n)= \frac{1}{n+1} \begin{pmatrix} 2n \\ n \end{pmatrix}$

电影购票

电影票每张50元，如果有$m+n$个人排队买票，其中$m$个人各持有50元面值的钞票1张，另外$n(n\le m)$个人各持有100元面值的钞票1张，而票房没有预备找零．有多少种方法可以将这$m+n$个人排成一列，顺序购票，使得无需因为等待找零而耽误时间？

分析：当手持100元的人买票时，前面至少有一个手持50元的人已经买过票了。把手持50元的人买票看做进栈, 用$+1$表示。把手持100元的人买票看做出栈, 用$-1$表示。从而转化为求解合法出栈方式的问题(当$m=n$时，栈最终为空；当$m>n$时，栈内留下$m-n$个元素。)

答案: $\begin{pmatrix} m+n \\ m \end{pmatrix}- \begin{pmatrix} m+n \\ m+1 \end{pmatrix}$

街区对角线问题

有一个$n\times n$的网格，需要从左下角走到右上角，每一步均为向上或向右。求不越过对角线的单调路径的个数。

分析：从左下角走到右上角，总的步数为$2n$步，其中$n$步向右，$n$步向上。不穿过对角线的走法需要向右走的步数大于等于向上走的步数。将向右看做入栈，用$+1$表示；将向上看做出栈，用$-1$表示。则可以转化为合法出栈方法个数的问题。

答案： $C(n)= \frac{1}{n+1} \begin{pmatrix} 2n \\ n \end{pmatrix}$ ^note2

参考

^note1. 错误修改: $\lim_{x\rightarrow 0^+} \frac{1-\sqrt{1-4x}}{2x} =\lim_{x\rightarrow 0^+} \frac{1}{\sqrt{1-4x}}= 1$ 和 $\lim_{x\rightarrow 0^+} \frac{1+\sqrt{1-4x}}{2x} = + \infty $. ↩

^note2. 不过对角线的走法，也可以从左上部分走。走法要求向上的次数大于等于向右的次数。条数与从对角线下方走一样，所以答案为$2C(n)$。参考http://www.cnblogs.com/wuyuegb2312/p/3016878.html ↩

矩阵连乘最小计算次数

发表于 2019-04-10 | 分类于 code

问题

$A_{2\times 10}, B_{10 \times 50}, C_{50 \times 20}$是三个矩阵，下标表示矩阵大小。现计算$A_{2\times 10} * B_{10 \times 50} * C_{50 \times 20}$。

$({\color{red}{AB}})_{2 \times 50}C_{50 \times 20}$的计算量是${\color{red}{2 *50 *10}}+2 *20 *50=3000$
$A_{2 \times 10}({\color{red}{BC}})_{10 \times 20}$的计算量是 ${\color{red}{10 *20 *50}}+2 *20 *10=10400$

可以观察到，不同的计算次序(在不同的地方加括号)有不同的计算量。

推广到一般情形：
给定$A=A_0 * A_1* … * A_{n-1}$，其中$A_i$的大小为: $d_i \times d_{i+1}$。选择计算顺序，使得(乘法)计算次数最小。

代码

def matrix_chain(d):
    ''' 
    d is a list: [d[0],d[1],...,d[n-1],d[n]]
    return: the total number of scalar of multiplications
    '''
    n = len(d)-1
    N = [[0]*n for _ in range(n)]
    for b in range(1,n):      # number of  products in subchain
        for i in range(n-b):  # start of subchain
            j = i+b
            N[i][j] = min([N[i][k]+N[k+1][j]+d[i]*d[k+1]*d[j+1] for k in range(i,j)])
    return N[0][n-1]

测试

1	matrix_chain([2,10,50,20])

Out: 3000

计算复杂度分析：$b,j,k$三重循环，所以为 $O(n^3)$

分析

使用动态规划(dynamic programming,DP)：确定迭代公式和确定边界条件

迭代公式

设$N_{0,i}$为$A=A_0 * … *A_i$的最小计算次数

计算$A_0 * A_1* … * A_{n-1}$一定会出现这种形式$(A_0 * … A_i) * (A_{i+1} * … * A_{n-1})$。

那么$(A_0 * … A_i) * (A_{i+1} * … * A_{n-1})$的最小计算次数为：

$\overbrace{N_{0,i}}^{(A_0 * ...* A_i)}+\overbrace{N_{i+1,n-1}}^{(A_{i+1} * ... * A_{n-1})}+\overbrace{d_0 * d_{i+1} * d_n}^{(A_0...A_i)*(A_{i+1}...A_{n-1})}$

$A$的最最小计算次数为$i$从0到n-2遍历后得到的最小计算量：

$N_{0,n-1} = \min_{i=0,1,...,n-2} \lbrace N_{0,i} + N_{i+1,n-1}+d_0 * d_{i+1} * d_n \rbrace$

其中

$N_{0,i}$的计算类似于$N_{0,n-1}$, 具体如下：
$N_{0,i} = \min_{j=0,1,...,i-2} \lbrace N_{0,j} + N_{j+1,i}+d_0 * d_{j+1} * d_{i+1} \rbrace$
$N_{i+1,n-1}$的计算公式如下：
$N_{i+1,n-1} = \min_{k=i+1,i+2,...,n-2} \lbrace N_{i+1,k} + N_{k+1,n-1}+d_{i+1} * d_{k+1} * d_n \rbrace$

这样就不断迭代到最基本的情况$N_{i,i}$(计算$A_i$) 和 $N_{i,i+1}$(计算$A_i A_{i+1}$)。

$A_i*A_{i+1} \cdots* A_{j}$的最优计算次数为
$N_{i,j} = \min_{k=i,i+1,...,j-1} \lbrace N_{i,k} + >N_{k+1,j}+d_i * d_{k+1} * d_{j+1} \rbrace$

边界条件

$N_{i,i}=0$：$A_i$不需要乘法计算
$N_{i,i+1}=d_i*d_{i+1}*d_{i+2}$: 计算$A_i A_{i+1}$

过程展示

初始化$N_{i,i}=0, \quad i=0,1,…,n-1$
计算$N_{i,i+1}\quad i=0,1,…,n-2 $
$\begin{align} N_{0,1} &= \min_{k=0} \{N_{0,0}+N_{1,1}+d_0d_1d_2\} \\ N_{1,2} &= \min_{k=1} \{N_{1,1}+N_{2,2}+d_1d_2d_3\}\\ \cdots \\ N_{n-2,n-1} &= \min_{k=n-2} \{N_{n-2,n-2}+N_{n-1,n-1}+d_{n-2} d_{n-1}d_{n}\} \end{align}$
计算$N_{i,i+2} \quad i=0,1,…,n-3$
$\begin{align} N_{0,2} &= \min_{k=0,1} \{N_{0,0}+N_{1,2}+d_0d_1d_3,N_{0,1}+N_{2,2}+d_0d_2d_3\} \\ \cdots \\ N_{n-3,n-1} &= \min_{k=n-3,n-2} \{N_{n-3,n-3}+N_{n-2,n-1}+d_{n-3} d_{n-2}d_{n}, \\ & \quad \qquad \quad \qquad N_{n-3,n-2}+N_{n-1,n-1}+d_{n-3} d_{n-1}d_{n}\} \end{align}$
$\cdots$
计算$N_{0, n}$

二叉树前序、中序和后序遍历

发表于 2019-01-20 | 分类于 code

二叉树(binary tree)有三种常见的遍历方式

前序遍历(preorder traversal): VLR

先访问树的根节点，然后使用前序遍历递归地访问根节点的左子树，最后使用前序遍历递归地访问根节点的右子树
中序遍历(inorder traversal): LVR

先使用中序遍历递归地访问根节点的左子树，然后访问根节点，最后使用中序遍历递归地访问根节点的右子树
后序遍历(postorder traversal): LRV

先使用后序遍历递归地访问根节点的左子树，然后使用后序遍历递归地访问根节点的左子树，最后访问根节点

这三种访问方式都是深度优先搜索(deep first search, DFS): 沿着某个方向一直访问节点，当到了根节点时(无法继续访问节点)，就返回到某个合适的节点继续沿着该固定方向访问。

下面二叉树的前序遍历VLR为：[0,1,3,4,2,5,6]

From http://mishadoff.com/blog/dfs-on-binary-tree-array/

本文用递归和非递归方式实现三种不同的遍历，完整代码见binary tree and its traverse。

构造树

# 定义节点
class node(object):
    def __init__(self,val=None):
        self.left = None
        self.right = None
        self.val = val
# 定义上图中的树
r = node(0)
r.left, r.right = node(1),node(2)
r.left.left, r.left.right = node(3),node(4)
r.right.left,r.right.right = node(5),node(6)

递归方式

将遍历方式定义为一个函数，然后根据不同遍历的定义写代码。

比如前序遍历VLR用preorderTraverse(root)表示。那么

$先访问树的根节点，\\ 然后使用前序遍历递归地访问根节点的左子树，\\ 最后使用前序遍历递归地访问根节点的右子树.$

可以表示为

def preorderTraverse(root):
    访问root的值
    preorderTraverse(root.left) #前序遍历左子树
    preorderTraverse(root.right) #前序遍历右子树

前序遍历VLR

def preorderTraverse(r):
    if r is None:
        return []
    result = []
    if r:
        result += [r.val]
    if r.left:
        result += preorderTraverse(r.left)
    if r.right:
        result += preorderTraverse(r.right)
    return result

输出

1	preorderTraverse(r)

Out:[3, 1, 4, 0, 5, 2, 6]

中序遍历LVR

def inorderTraverse(r):
    if r is None:
        return []
    result = []
    if r.left:
        result += inorderTraverse(r.left)
    if r:
        result += [r.val]
    if r.right:
        result += inorderTraverse(r.right)
    return result

输出

1	inorderTraverse(r)

Out:[3, 1, 4, 0, 5, 2, 6]

后序遍历LRV

def postorderTraverse(r):
    if r is None:
        return []
    result = []
    if r.left:
        result += postorderTraverse(r.left)
    if r.right:
        result += postorderTraverse(r.right)
    if r:
        result += [r.val]
    return result

输出

1	postorderTraverse(r)

Out:[3, 4, 1, 5, 6, 2, 0]

非递归方式

非递归方法比较复杂一些。

深度优先搜索(DFS)算法是一种使用回溯思想的递归算法: 它首先沿着某方向一直搜索直到无可访问节点，然后通过回溯继续访问其它节点。

这类过程适合用数据结构中的栈stack来实现。分为如下3个过程：

选择一个起始节点并将其所有相邻节点推入堆栈。
从堆栈弹出一个节点作为访问的下一个节点，并将其所有相邻节点推入堆栈。
重复此过程，直到堆栈为空。注意，需要标记或者通过逻辑防止二次访问某节点。如果程序处于无限循环中，可检查是否有标记或者逻辑错误。

前序遍历VLR

def preorderTraverse2(r):
    if r is None:
        return []
    stack = []
    result = []
    node = r   # 起始节点
    while stack or node:
        # 一直访问左子树直到根节点，同时将访问过的节点添加到stack中,
        # 以便回溯时使用
        while node:
            result.append(node.val)
            stack.append(node)
            node = node.left
        # 回溯
        temp = stack.pop()
        node = temp.right
    return result

输出

1	preorderTraverse2(r)

Out:[3, 1, 4, 0, 5, 2, 6]

中序遍历LVR

def inorderTraverse2(r):
    if r is None:
        return []
    stack = []
    result = []
    node = r
    while stack or node:
        # 先走到最左的根节点
        while node:
            stack.append(node)
            node = node.left
        # 然后开始访问节点的值
        temp = stack.pop()
        result.append(temp.val)
        node = temp.right
    return result

输出

1	inorderTraverse2(r)

Out:[3, 1, 4, 0, 5, 2, 6]

后序遍历LRV

将后序遍历LRV,看做VRL的逆。
所以，先对做VRL遍历，然后翻转节点列表，得到后序遍历列表。

def postorderTraverse2(r):
    if r is None:
        return []
    stack = []
    result = []
    node = r
    # VRL
    while stack or node:
        while node:
            result.append(node.val)
            stack.append(node)
            node = node.right
        temp = stack.pop()
        node = temp.left
    # reverse VRL to get LRV
    result = result[::-1]
    return result

输出

1	postorderTraverse2(r)

Out:[3, 4, 1, 5, 6, 2, 0]

参考

Problem Solving with Algorithms and Data Structures-6.7. Tree Traversals
DFS on Binary Tree Array
Depth First Search

字符串编辑距离-edit distance of two strings

发表于 2019-01-18 | 分类于 code

问题

给2个字符串str1和str2. 通过对单个字符进行插入、删除或替换的方式将str1变换为str2的最小操作次数，定义为编辑距离。

输入任意字符串str1和str2，输出将str1变换为str2的编辑距离。

LeetCode：72. Edit Distance

Given two words word1 and word2, find the minimum number of operations required to convert word1 to word2. You have the following 3 operations permitted on a word:

Insert a character

Delete a character

Replace a character

Jupyter notebook: edit distance of two strings

代码

def editDistance(str1,str2):
    #异常处理
    if str1 is None or str2 is None:
        return -1 
    m,n = len(str1),len(str2)
    # 定义编辑距离
    d = [[0]*(n+1) for _ in range(m+1)]
    # 边界条件
    for i in range(1,m):
        d[i][0] = i
    for j in range(1,n):
        d[0][j] = j
    # 迭代过程    
    for i in range(1,m+1):
        for j in range(1,n+1):
            if str1[i-1]==str2[j-1]:
                d[i][j] = d[i-1][j-1]
            else:
                d[i][j] = min(d[i][j-1]+1,    # 插入
                              d[i-1][j]+1,    # 删除
                              d[i-1][j-1]+1)  # 替换
    return d[-1][-1]

测试

1
2
3

str1 ='cat'
str2 = 'mouse'
editDistance(str1,str2)

Out: 3

分析

使用动态规划(dynamic programming,DP)，一般DP有两个过程：确定迭代公式和确定边界条件

确定迭代公式

d[i][j]表示str1[0:i]与str2[0:j]之间的编辑距离，

str1[0:i]不包含指标(index)为i的字符。如果str1='insert', 则str1[0:3]='ins'

假设已计算出$\{d[p][q]: 0 \le p \le (i-1),0 \le q \le (j-1) \}$。有如下两种情况：

当前字符相等：str1[i-1]==str2[j-1]

此时，只需要将字符串str1[0:i-1]转换成str2[0:j-1], 即

d[i][j]=d[i-1][j-1]

比如 $ mouse_{i-1} t\rightarrow mouse_{j-1} \Rightarrow mous_{i-2} \rightarrow mous_{j-2} $

当前字符不等：str1[i-1]!=str2[j-1]

此时，就需要使用插入、删除或替换的方式处理

比如问题是 $ mouse_{i-1} \rightarrow cat_{j-1}$

插入: 在$e_{i-1}$后插入一个$t$, 问题转化为 $ mouse_{i-1} t\rightarrow cat_{j-1} \Rightarrow mouse_{i-1} \rightarrow ca_{j-2} $

从而：d[i][j]=d[i][j-1]+1

插入1次，$ mouse_{i-1} $修改为$ca_{j-2}$修改d[i][j-1]次

删除：删除$e_{i-1}$，问题转化为$ mous_{i-2} \rightarrow cat_{j-1}$

从而：d[i][j]=d[i-1][j]+1

删除1次，$ mous_{i-2} $修改为$cat_{j-1}$修改d[i-1][j]次

替换：把 $e_{i-1}$替换为$t$, 问题转化为 $ moust_{i-1} \rightarrow cat_{j-1} \Rightarrow mous_{i-2} \rightarrow ca_{j-2} $

从而: d[i][j]=d[i-1][j-1]+1

替换1次，$ mous_{i-2} $修改为$ca_{j-2}$修改d[i-1][j-1]次

第2种情况下，d[i][j]取d[i-1][j]+1,d[i][j-1]+1,d[i-1][j-1]+1中最小值。

d[i][j] = min(d[i-1][j]+1,d[i-1][j-1]+1,d[i][j-1]+1)

边界条件

当str1=' '时，使用插入将str1变换为str2, d[0][j]=j

' '变为 cat, 编辑次数为3次

当str2=' '时，使用删除将str1变换为str2, d[i][0]=i

mouse变为' ',编辑次数5次

批标准化-batch normalization

发表于 2019-01-08 | 分类于 deep learning

上一篇标准化-normalization中说到对线性回归模型，当输入特征在各个维度上取值范围相差较大时，可通过对输入标准化(normalization)来加快训练速度。

其中用到的回归模型$y = w_1x_1+w_2 x_2$可以看做bias为0，激活函数为线性函数的单个神经元模型。

借鉴输入标准化可以加快训练速度的思想，批标准化(Batch-normalization, BN) 通过对深度网络每一层每个连接单元标准化来加快网络的训练速度。

Batch-normalization(BN)来自于文章Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift .

本文希望回答下面几个问题：

(1) 批标准化(BN)如何做的？

(2) BN除加快训练速度外，还有什么优点？

模型

一个多层神经网络模型：

其中

$x$是输入
$W^1,W^2,\cdots$是权重
$z$ 为激活函数的输入
$a$ 为激活函数的输出

标准化(normalization) 就是对每层激活函数的输入做标准化。设整个数据集的均值和方差分别为$\mu_{}$和$\sigma_{}$，那么标准化公式为：

$\tilde{z} = \frac{z-\mu_{}}{\sigma_{}}.$

批标准化(batch-normalization)

但是多数深度学习模型数据量比较大，内存不足以计算整个数据集的均值和方法；或者为在线学习(online learning)，数据动态获取，无法得到整个数据集。

对这种情况，就需要用batch-normalization：使用批量数据(样本)的均值和方差替代整体(总体)的均值和方差。

既然是用样本估计总体，batch选的越大，估计值就越准。所以batch选的越大越好。

具体就是使用batch数据的得到的均值和标准差的移动平均值，近似整个数据的均值和标准差。

另外为了增加网络的表示能力^表示能力，对标准化后的数据加了偏移项$\beta$及尺度参数$\gamma$ ($\beta,\gamma$的维数与$\tilde{z}$一样):

$\hat{z} = \gamma \odot \tilde{z}+\beta$

其中$\odot$是分量与分量的乘积(element-wise product)。

问题(1)的回答

将数据集$D$分为训练集$D^{train}=(X^{train},y^{train})$和测试集$D^{test}=(X^{test},y^{test})$.

训练阶段:

从训练集$D^{train}=(X^{train},y^{train})$中从抽样一个batch-size的数据$\{x^i,y^i\}_{i=1}^{bs}$, 然后在每个激活函数前做标准化+平移+变化尺度：

image modified based on Machine Learning and having it deep and structured (2018,Spring): Special Training Technology

其中

$u_{bs},\sigma_{bs}$与每个batch的数据有关
$\beta,\gamma$作为模型的参数通过误差反向传播求解
第一个batch数据用$u_{bs},\sigma_{bs}$，后面每个batch数据用下面描述的$\bar{\mu}^{t}$和$\bar{\sigma}^{t}$。

在迭代过程中记录其指数移动平均(exponential moving average)值$\bar{\mu}^{t}$和$\bar{\sigma}^{t}$，用来近似整个数据的均值$\mu$和方差$\sigma$。

设共迭代$T$步:

当$t=1$时，

$\begin{align} \bar{\mu}^{1}=\mu_{bs}\\ \bar{\sigma}^{1}=\sigma_{bs} \end{align}$

当$t\ge2$时,

$\begin{align} \bar{\mu}^{t+1} &= \alpha_1 \mu^t_{bs} +(1-\alpha_1) \bar{\mu}^t_{} \tag{1} \\ \bar{\sigma}^{t+1} &= \alpha_2 \sigma_{bs}^t +(1-\alpha_2)\bar{\sigma}_{}^t \tag{2} \end{align}$

其中

$\mu^t_{bs}$-第$t$步迭代一个Batch的平均值；
$\bar{\mu}^t_{}$-第$t$步迭代模型中使用的均值；
$0 < \alpha_1< 1$-权重延迟系数。$\alpha_1$越大，过去均值对现在均值影响越小；
公式(2)与公式(1)类似。

将最后一步迭代的参数值记为$\mu_{train}=\bar{\mu}^{T},\sigma_{train}=\bar{\sigma}^{T}$.

测试阶段:

计算测试集输入$X^{test}$的均值和标准差分别为$\mu_{test}$, $\sigma_{test}$，则测试过程的均值和方差$\mu,\sigma$取值如下:

$\begin{align} \mu &= \alpha_1 \mu_{test} +(1-\alpha_1) \mu_{train} \\ \sigma &= \alpha_2 \sigma_{test} +(1-\alpha_2)\sigma_{train} \end{align}$

Internal Covariate Shift

文章Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift 将各层参数分布不一样而导致激活值(activations)分布不一样的现象定义为Internal Covariate Shift ^{Internal Covariate Shift}.

为观察Internal Covariate Shift 现象，对cifar10分类问题，使用keras搭建一个CNN模型：12层的网络, 其中有9个卷积层。代码见jupyter notebook。

因为activations直接与参数有关，所以我们通过可视化参数的分布来说明activation的分布变化。

层数太少，比如5层，并没有观测到分布的变化现象。

训练结束后，权重的分布

没有批标准化(左图), conv1 conv2 conv3 差别还是较大的，其它卷积层权重分布差不多
有批标准化(右图), conv1 与其他卷积层差别较大，其它层权重分布相对比较一致

其它优点

Batch-normalization最大的优点是加快训练速度，除此之外，它还有一些其它的优点:

问题(2)的回答

(1) 避免梯度消失(gradient vanish)

网络中使用 sigmoid 或者 tanh 激活函数时，批标准化可调整输入到0附近，避免了数据落入导数为0的区域，从而可以一定程度上避免梯度消失问题。

(2) 正则化(regularization)

给定网络的损失函数: $\min_w L(w)=\sum (y_i-\hat{y}_i(w))^2+|w|_2$

使用批标准化后, $w’=kw(k \geq 1)$ 与 $w$的预测值是一样的(如下图)

但是 $L(w’) \geq L(w)$, 所以训练中会选择较小的参数值$w$，从而达到正则化的效果。

^{Internal Covariate Shift}. We define Internal Covariate Shift as the change in the distribution of network activations due to the change in network parameters during training. 来自文章BN:… by Reducing Internal Covariate Shift ↩

^表示能力. 标准化将网络输入值变换到一定的范围内，限制了网络的表示能力。比如激活函数为Relu时，大部分输入将被限制在$-$3到$3$以内，其它值出现的概率较小。 ↩

使用tensorflow,opencv,scipy实现高斯模糊

发表于 2018-12-29 | 分类于 machine learning

概念

图像的高斯模糊就是图像与高斯权重(高斯核)做卷积。

简单来说，就是把每个像素点变成周围像素点的加权平均。

对比

本文展示三种不同软件包tensorflow,opencv,scipy实现高斯模糊的方法。代码见Gaussian blur-tensorflow,opencv,scipy.

tensorflow版本可以嵌入到深度学习模型的代码中，方便求导等运算。

opencv 相对scipy方便的多：灰度图和彩色图可使用同一函数，并且可直接指定高斯核大小。

几种实现方法的简要对比:

	灰度图（单通道)	彩色图(多通道)	备注
tensorflow	tf.nn.conv2d^tf	tf.nn.conv2d	多通道图分别卷积，然后合并
opencv	cv.GaussianBlur^cv	cv.GaussianBlur
scipy	filters.gaussian_filter^filters	filters.gaussian_filter	多通道图分别卷积，然后合并。核大小只能间接指定。

过程

高斯核的计算分成两步: (1) 计算正态权重 (2) 归一化正态权重得到高斯核。

正态权重$G$的计算公式

$G(x,y)=\frac{1}{2\pi\sigma^2}e^{-\frac{x^2+y^2}{2\sigma^2}}$

其中

$\sigma$是模糊因子(blurring factor). 越大模糊的越厉害；
$x$是水平方向上像素点离中心点的距离；
$y$是竖直方向上像素点离中心点的距离。

一个$3\times3$的卷积核，假设其中心点坐标为$(x_c,y_c)=(0,0)$, 则其相邻像素点坐标$(x,y)$如下:

$\begin{bmatrix} (-1,1) & (0,1) & (1,1) \\ (-1,0) & \color{red}{(0,0)} & (1,0) \\ (-1,-1) & (0,-1) & (1,-1) \end{bmatrix}$

假设$\sigma=1.5$, 则$G$的计算结果如下:

gaus

image from How to program a Gaussian Blur without using 3rd party libraries

对$G$做归一化，可得高斯核$K$为(保留2位小数):

$\begin{align} K &= G/sum(G)=G/0.48 \\ &= \begin{bmatrix} 0.09 & 0.12 & 0.09 \\ 0.12 & \color{red}{0.15} & 0.12 \\ 0.09 & 0.12 & 0.09 \end{bmatrix} \end{align}$

通用高斯核的python实现

def gaussian_kernel(size=5,sigma=2):
    '''
    size: int,一般取为奇数
    sigma: blur factor
    
    return: (normalized) Gaussian kernel，大小 size*size
    '''
    x_points = np.arange(-(size-1)//2,(size-1)//2+1,1)
    y_points = x_points[::-1]
    xs,ys = np.meshgrid(x_points,y_points)
    kernel = np.exp(-(xs**2+ys**2)/(2*sigma**2))/(2*np.pi*sigma**2)
    return kernel/kernel.sum()

可视化高斯核

1 2	plt.imshow(gaussian_kernel(20,5),cmap='Reds') plt.colorbar()

灰度图高斯模糊

将图形与$K$做卷积，得到模糊图。

彩色图(RGB)高斯模糊

将原图的3个通道分别与高斯核$K$做卷积，得到不同通道的模糊图，然后合并得到最终模糊图。

实现

tensorflow

灰度图

计算kernel的值，然后将其扩展为tensorflow可用的维度

1 2	kernel = gaussian_kernel(size=10,sigma=5) kernel = kernel[:, :, np.newaxis, np.newaxis] # height,width, channel_in, channel_out

灰度图高斯卷积

1 2	x = tf.placeholder(tf.float32,shape=(None,200,300,1)) # (batch_size, height, width, channel) x_blur = tf.nn.conv2d(x, kernel, strides=[1, 1, 1, 1], padding='SAME')

计算模糊图(图像image_gray也需要调整为tensorflow可用的维度

s = tf.InteractiveSession()
image_batch = image_gray[np.newaxis,:,:,np.newaxis]  # (batch_size, height, width, channel)
blur_imgs = s.run(x_blur,feed_dict={x:image_batch})
s.close()

可视化结果

彩色图

第2步需要将原图分开为单通道图xr,xg,xb, 然后分别做卷积得到xr_blur,xg_blur,xb_blur, 最后合并得到彩色图的模糊图xrgb_blur:

x = tf.placeholder(tf.float32,shape=(None,200,300,3))    # (batch_size, height, width, channel)
xr,xg,xb =tf.expand_dims(x[:,:,:,0],-1),tf.expand_dims(x[:,:,:,1],-1),tf.expand_dims(x[:,:,:,2],-1)

xr_blur = tf.nn.conv2d(xr, kernel, strides=[1, 1, 1, 1], padding='SAME')
xg_blur = tf.nn.conv2d(xg, kernel, strides=[1, 1, 1, 1], padding='SAME')
xb_blur = tf.nn.conv2d(xb, kernel, strides=[1, 1, 1, 1], padding='SAME')
xrgb_blur = tf.concat([xr_blur,xg_blur,xb_blur],axis=3)

结果:

opencv

使用cv.GaussianBlur。灰度图和彩色图都使用该函数。

可以指定高斯核大小及模糊因子大小。

import cv2 as cv

size = 5
sigma = 5
blur = cv.GaussianBlur(image,(size,size),sigma)

scipy

使用scipy.ndimage.filters.gaussian_filter.

灰度图可直接使用gaussian_filter, 彩色图需要分别对不同通道做模糊处理。

该函数通过truncate参数指定高斯核的大小$s$，truncate与$s$和模糊因子$\sigma$的关系是truncate$=(s-3)/ (2\sigma)$^解释.

高斯模糊函数定义为

def gaussian_blur(image, sigma=2, size = 3):
    '''
    image: [h,w,c]
    size: 高斯核大小
    sigma：模糊因子
    '''
    # compute the truncate using size
    t = (((size - 1)/2)-0.5)/sigma
    im = np.zeros_like(image)
    channels = image.shape[-1]
    for i in range(channels):
        im[:, :, i] = filters.gaussian_filter(image[:, :, i], sigma,truncate=t)
    return im

结果:

^解释. 参考https://stackoverflow.com/questions/25216382/gaussian-filter-in-scipy ↩

^tf. import tensorflow as tf ↩

^cv. import cv2 as cv ↩

^filters. from scipy.ndimage import filters ↩

标准化-normalization

发表于 2018-12-18 | 分类于 machine learning

本文希望回答下面的问题：

对输入数据标准化为什么能加快训练速度？

所有过程代码见 juputer notebook.

线性回归

给定数据集$\{x^i,y^i\}_{i=1}^N$, 其中$x^i=(x^i_1,x^i_2) \in \mathbb{R}^2$,$y^i \in \mathbb{R}$。

回归模型 $y = w^Tx=w_1x_1+w_2 x_2$ (设模型中的b为0）

损失函数取为均方误差

$L(w_1,w_2) = \frac{1}{N} \sum_{i=1}^N L_i(x^i,y^i) = \frac{1}{N} \sum_{i=1}^N (w_1x_1^i+w_2 x_2^i-y^i)^2 \tag{1}$

为了简洁，$L$也表示$L(w_1,w_2) $

损失函数在$w_1 = c$的截面记为 $L(w_1=c,w_2) $或者$L(w_2)$

损失函数在$w_2 = c$的截面记为 $L(w_1,w_2=c) $或者$L(w_1) $

数据模拟

数据点$\{x^i,y^i\}_{i=1}^{200}$由线性回归 $y=0.1 x_1+0.2 x_2+\epsilon$模拟得到。

其中（$x_2$的取值范围是$x_1$的10倍)

$x_1 \sim 1*N(0,1)$
$x_2 \sim 10*N(0,1)$
$\epsilon$为噪音, 分布为 $\epsilon \sim N(0,1)$

损失函数

将200个数据点代入公式(1), 得到$L, L(w_1,w_2=0.1), L(w_1=0.2,w_2)$的几何图形:

左图：损失函数在参数平面($w_1,w_2$)上具有椭圆形水平集。

中图：损失函数沿着$w_2$的截面图

右图：损失函数沿着$w_1$的截面图

可以观察到损失函数$L$在不同参数方向上的下降速度(梯度)差别很大，这会导致什么问题呢？

在$w_1$的方向上选择最佳学习率${\eta_{w_1}}$更新参数 $w_1\leftarrow w_1-\eta_{w_1} \nabla_{w_1} L$，

在$w_2$的方向上使用$\eta_{w_1}$更新参数 $w_2\leftarrow w_2-\eta_{w_1} \nabla_{w_2} L$.

结果是在$w_2$方向上由于学习率太大，会发散。
在$w_2$的方向上选择最佳学习率$\eta_{w_2}$更新参数 $w_2\leftarrow w_2-\eta_{w_2} \nabla_{w_2} L$，

在$w_1$的方向上使用$\eta_{w_2}$更新参数$w_1\leftarrow w_1-\eta_{w_2} \nabla_{w_1} L$.

结果是在$w_1$方向上由于学习率太小，收敛速度太慢。

所以，为了避免在其中一个方向上发散，只能选择较小的学习率$\eta_{w_2}$.

$L$在$w_1$和$w_2$下降速度不一样，具体表现在$\nabla_{w_2} L$比$\nabla_{w_1} L$大很多，从而在整个损失函数上产生如下锯齿形的更新过程, 从而收敛速度较慢。

下图是梯度下降算法最小化损失函数的参数更新过程：

标准化

为了解决这个问题，加快收敛速度：

一种办法是自动调整不同方向的学习率，在$w_1$方向上把学习率调大，在$w_2$方向上把学习率调小。从而诞生了许多自适应学习率算法Adagrad, RMSprop, Adam等等。
另一种解决方法就是调整输入的尺度(scaling), 使其归一化到同一个范围^{featureScaling}，比如[-1,1]; 或者是通过标准化^{normalization}，转化为0均值，单位方差的输入。

采用2的方法。对输入特征$\{x^i\}_{i=1}^N$做标准化$\tilde{x}^i = \frac{x^i-\mu}{\sigma}$ ,其中

$\mu = [\mu_1,\mu_2]^T $, $\mu_j = \frac{1}{N} \sum_{i=1}^N x_j^i\quad j=1,2$
$\sigma = [\sigma_1, \sigma_2]^T$, $\sigma_j= \frac{1}{N} \sum_{i=1}^N (x^i_j-\mu_j)^2 \quad j=1,2$

然后重新计算损失函数，可得下图：

下图是梯度下降算法最小化损失函数(输入标准化)的参数更新过程：：

观察图形，可以得到下面重要的结论：

使用标准化(normalization)后，损失函数从椭圆形改进到接近圆形；
$L$在$w_1$和$w_2$下降速度差的不多, 所以$\nabla_{w_2} L$和$\nabla_{w_1} L$也差不多, 在$w_1,w_2$方向上的学习率也差不多, 从而可以使用较大的学习率。这样就避免了锯齿形的更新方式, 加快了收敛速度。

下图是训练误差图：

(1). 没有标准化，学习率lr为0.011时收敛，0.0111时发散

(2). 标准化后，学习率lr可以取得更大(比如，为未标准化的10倍左右)

(3). 标准化后迭代到20步基本收敛到最优值，未标准化迭代到100步仍未收敛到最优

^{featureScaling}. $\frac{x-x_{min}}{x_{max}-x_{min}}$ ↩

^{normalization}. $\frac{x-\mu}{\sigma}$ ↩

主成分分析-principal component analysis

发表于 2018-12-10 | 分类于 machine learning

问题描述

某些高维空间的数据可以嵌入或者近似嵌入在低维空间中。

例子：$\mathbb{R}^2$空间的数据点$\{(1,1),(2,2),…,(5,5)\}$可以看做$\mathbb{R}^1$空间的点$\{1,2,..,5\}$

我们的目的是将一些数据点$\{x_1,x_2,…,x_N\} , x_i \in \mathbb{R}^p$ 变换成$\{y_1,y_2,…,y_N\} , x_i \in \mathbb{R}^q$，其中$q \le p$, 另外要求$\{y_1,y_2,…,y_N\}$的样本协方差矩阵是对角阵。

随机变量才有方差或者协方差，在下面的推导中，我们

将$\{x_1, x_2, …, x_N\}, x_i \in \mathbb{R}^p$ 理解为来自随机变量$x$的样本
将$\{y_1, y_2, …, y_N\}, y_i \in \mathbb{R}^q$ 理解为来自随机变量$y$的样本

$x$的(总体)均值和(总体)协方差分别用$Ex$和$Var(x)=cov(x,x)$表示。

$x$的样本均值和样本协方差分别用$\bar x =\frac{1}{N}\sum_{n=1}^N x_n$和$\frac{1}{N}\sum_{n=1}^N(x_n-\bar x)(x_n-\bar x)^T$表示。

问题分析

高维空间的点，可以通过投影得到低维空间的点。

考虑投影到一维空间$\{ ku_1|k \in \mathbb{R}，u_1 \in \mathbb{R}^p \}$，我们希望投影数据能包含原始数据更多的信息。

信息的度量就等于不确定性的多少, 具体可以用熵(entropy)来量化, 熵越大表示信息越多^智能时代。

假设投影数据服从正太分布 $N(\mu,\sigma^2)$, 其熵(entropy)为 $\frac{1}{2} \ln(2\pi e \sigma^2)$.

从而，分布的方差$\sigma^2$越大，熵也就越大，包含的信息也就越多。

所以，我们需要选择投影数据方差最大的方向$u_1$。

另一种理解方式是最小化投影误差(projection error minimization)^PRML.

计算过程

先分析最简单的情形: 把数据投影到一维空间，然后将一维情形得到的结果推广到高维空间。

投影到一维: $\mathbb{R}^p \Rightarrow\mathbb{R}^1$

因为投影数据的方差只与$u_1$的方向有关，而与其大小无关。

所以为了便于计算，我们限制$u_1$为单位向量，即$u_1^T u_1=1$。

计算$\{x_1,x_2,…,x_N\} $投影在$u_1$上的投影$\{y_1,y_2,…,y_N\} $。
$x_n$在$u_1$上投影为:
$y_n = \frac{u_1^T x_n}{u_1^T u_1} = u_1^T x_n$
计算投影数据 $\{y_1,y_2,…,y_N\}$的方差：
$\begin{align} Var(y) &= E(y-Ey)^2 \\ &\approx \frac{1}{N}\sum_{n=1}^N(u_1^Tx_n - u_1^T \bar{x})^2 \leftarrow \bar{x}=\frac{1}{N} \sum_{n=1}^N x_n \\ &=\frac{1}{N}\sum_{n=1}^N (u_1^Tx_n - u_1^T \bar{x})(u_1^Tx_n - u_1^T \bar{x}) \\ &= \frac{1}{N}\sum_{n=1}^N (u_1^Tx_n - u_1^T \bar{x})(u_1^Tx_n - u_1^T \bar{x})^T \leftarrow a=a^T \text{ if a is a scaler} \\ &= u_1^T \left[ \frac{1}{N}\sum_{n=1}^N (x_n - \bar{x})(x_n - \bar{x})^T \right]u_1 \leftarrow S=\frac{1}{N}\sum_{n=1}^N (x_n - \bar{x})(x_n - \bar{x})^T \\ &=u_1^T S u_1 \end{align}$
上式中$S$是$\{x_1,x_2,…,x_N\} $的样本协方差矩阵。
最大化方差，求解$u_1$

目标函数：
$\max_{u_1} \qquad u_1^T S u_1 \\ \quad \text{s.t.} \qquad u_1^T u_1=1$
这是一个有限制条件的优化问题，可以通过拉格朗日乘子法将其转化为无约束优化问题：
$\max_{u_1} \quad F(u_1)=u_1^T S u_1 +\lambda_1 (1- u_1^T u_1).$F(u_1)$取到最大值满足的条件：$
$F(u_1)$取到最大值满足的条件：
$\begin{align} \frac{\partial F(u_1) }{\partial u_1} &= 2Su_1-2\lambda_1 u_1=0 \\ 1- u_1^T u_1&=0 \end{align}$
从而，解$u_1$就是$S$的最大特征值$\lambda_1$对应的特征向量。此时，$F(u_1)=\lambda_1$.

$u_1$称为第一主成分(first principal component), 记为pc1。得到的一维空间 $\mathbb{R}^1 = \{ ku_1|k \in \mathbb{R}，u_1 \in \mathbb{R}^p \}$称为第一主子空间(first principal subspace)。

结果:

投影矩阵为$u_1^T$，投影数据为 $y_i = u_1^Tx_i \in \mathbb{R}^1 \qquad i =1,2,…,N$;
$y = u_1^Tx \in \mathbb{R}^1$，其方差为 $\lambda_1$。

投影到多维: $\mathbb{R}^p\Rightarrow\mathbb{R}^q(2\le q\le p)$

当$p=1$时, 我们已经说明了保持方差最大的投影方向是$S$的最大特征值$\lambda_1$对应特征向量$u_1$的方向。

当$p=M$时, 假设保持方差最大的投影方向依次是$S$的前$M$个特征值$\lambda_1 \ge \lambda_2 \ge …\ge\lambda_M$对应的特征向量$u_1,u_2,…,u_M$的方向。

$u_1,u_2,…,u_M$分别是第一主成分(first principal component), …,第M个主成分(M-th principal component). 记为pc1,pc2,…,pcM.

$u_1,u_2,…,u_M$构成的空间$\mathbb{R}^M = \{ k_1u_1+…+k_M u_M|k_i \in \mathbb{R}，u_i \in \mathbb{R}^p\}$称为M维主子空间(M dimensional principal subspace)。

现证明$p=M+1$时，结论仍然成立。

此时，优化的目标函数为：

$\max_{u_{M+1}} \qquad u_{M+1}^T S u_{M+1}$

约束条件($u_{M+1}^T$为单位向量且与前$M$个主成分正交)为:

$\begin{align} u_{M+1}^T u_{M+1}&=1 \quad i=1,2,...,M\\ u_{M+1}^T u_i&=0 \quad i=1,2,...,M \end{align}$

同样，这是一个有限制条件的优化问题，可以通过拉格朗日乘子法将其转化为无约束优化问题：

$\max_{u_{M+1}} \quad F(u_{M+1})=u_{M+1}^T S u_{M+1}+\lambda_{M+1} (1- u_{M+1}^T u_{M+1})+\sum_{i=1}^M \eta_i u_{M+1}^T u_i.$

$F(u_{M+1})$取到最大值满足的条件：

$\begin{align} \frac{\partial F(u_{M+1}) }{\partial u_{M+1}} &= 2Su_{M+1}-2\lambda u_{M+1}+\sum_{i=1}^M \eta_i u_i =0 \tag{1}\\ 1- u_{M+1}^T u_{M+1}&=0 \\ u_{M+1}^T u_{i}&=0 \quad \text{i=1,2,...,M} \end{align}$

在公式(1)上乘以$u_{j}^T, j \in \{1,2,…,M\}$, 得到:

$u_{j}^Tu_{M+1}=0$.

$j=i$时，$ u_{j}^Tu_i = 1 $；否则$ u_{j}^Tu_i = 0$.

$S$是样本协方差矩阵，是对称的。$S^T=S$.

$u_{j}^TSu_{M+1}$ 是一个标量，所以$u_{j}^TSu_{M+1}=(u_{j}^TSu_{M+1})^T=u_{M+1}^TSu_j =\lambda_i u_{M+1}^T u_j$.

$\begin{align} 0&=2u_{j}^TSu_{M+1}-2\lambda u_{j}^Tu_{M+1}+\sum_{i=1}^M \eta_i u_{j}^Tu_i \\ &=2u_{M+1}^TSu_j-0+\eta_j \leftarrow Su_j=\lambda_j u_j \\ &=2\lambda_j u_{M+1}^Tu_j-0+\eta_j \\ &= 0+\eta_j\\ &=\eta_j \end{align}$

将$\eta_j=0$代入公式(1)，得到$Su_{M+1}=\lambda u_{M+1}$.

因为前$M$个最大的特征值已经使用过，所以只能选择第$M+1$大的特征值$\lambda_{M+1}$及其对于的特征向量$u_{M+1}$.

从而第$M+1$个主成分为$u_{M+1}$。

$M+1$维主子空间($M+1$ dimensional principal subspace) 为 $\mathbb{R}^{M+1} = \{ k_1u_1+...+k_M u_M+k_{M+1}u_{M+1}|k_i \in \mathbb{R}，u_i \in \mathbb{R}^p\}$ .

上面使用数学归纳法证明了数据$\{x_1,x_2,…,x_N\} , x_i \in \mathbb{R}^p$的$q$维主子空间($q$ dimensional principal subspace), 由数据的样本协方差矩阵$S=\frac{1}{N}\sum_{n=1}^N (x_n - \bar{x})(x_n - \bar{x})^T$的前$q$个依次减小的最大特征值$\lambda_1 \ge … \ge \lambda_q$对应的特征向量$u_1,…,u_q$构成。

结果：

变换矩阵为$U^T$, 投影数据为:

$y_i = \left[ \begin{matrix} u_1^T x_i \\ \cdots \\ u_{q}^T x_i \end{matrix} \right] = U^Tx_i \in \mathbb{R}^q \qquad i=1,2,...,N;$

$y=U^Tx$, $y$的协方差矩阵为:

$\begin{align} Var(y)&=cov(y,y) = cov(U^Tx,U^Tx) \\ &=U^Tcov(x,x)U = U^T S U\\ &= \left[ \begin{matrix} u_1^T \\ \cdots \\ u_{q}^T \end{matrix} \right] S \left[ \begin{matrix} u_1 & \cdots &u_{q} \end{matrix} \right]\\ &= \left[ \begin{matrix} u_1^TSu_1 & \cdots & u_1^T S u_q \\ & \cdots & \\ u_q^T S u_1 & \cdots & u_q^T S u_q \end{matrix} \right] \\ &= \left[ \begin{matrix} \lambda_1 & 0 & 0 \\ & \cdots & \\ 0 & 0 & \lambda_q \end{matrix} \right]. \end{align}$

一句话描述

对随机变量$x \in \mathbb{R}^p $, 寻找一个线性变换$y_{q \times 1}=W_{q\times p}x_{p \times 1} \in \mathbb{R}^q$, 使得变换后随机变量$y$的协方差矩阵是对角阵。

具体：设$\{x_1,x_2,…,x_N\} $为来自随机变量总体$x$的样本，其样本协方差矩阵为$S=\frac{1}{N}\sum_{n=1}^N (x_n - \bar{x})(x_n - \bar{x})^T$

$S$的特征值和特征向量关系：$S_{p\times p}U_{p \times q}=U_{p \times q}\Lambda_{q \times q}$

其中

$\Lambda$的对角元素为$S$的特征值，$\Lambda=\mathrm{diag}[\lambda_1..,\lambda_p]，\lambda_1\ge \lambda_2\ge… \lambda_p $；
$U$的相应列对应于$S$的特征向量, 分别称为第一、第二…第$p$个主成分。

此时，$W=U^T$. 即， $y=U^Tx$.

可视化

下面两个例子主要用来展示PCA降维的作用。

CS229: Machine Learning. Lecture 15-Principal Component Analysis.

分析了PCA的几个主要作用：

压缩(compression)-使用投影后的低维数据替代原始数据，用来聚类等；

降维(dimension reduce)-使用投影后的低维数据作为机器学习模型(线性回归，逻辑回归)的输入, 减少模型参数的个数，从而降低模型的复杂度，降低模型的方差(variance), 达到避免过拟合(overfitting)的目的；

去噪(noise reduction)-大特征值对应的主成分子空间表示主要的特征，而小特征值对应的主成分子空间表示次要的特征(在某些情形下, 就是噪音)。所以保留高特征值对应的主子空间的特征就达到了去除噪音的目的。

三维空间数据

左：原始数据 $\qquad \qquad \quad \quad$ 中：原始数据在第一主成分pc1和第二主成分pc2的投影

右上：原始数据取值范围 $\qquad $ 右下：投影数据取值范围

pca3d

image comes from http://setosa.io/ev/principal-component-analysis/

可从得到的结论：

PCA可以理解为坐标系的变换;
原始数据在各个主成分投影的方差：pc1最大，pc2次之，pc3最小。

MNIST数据集

MNIST是$28\times28$的图片，转化为向量后，每一张图片是空间$\mathbb{R}^{784}$上的点。

对该数据集做PCA，将其投影到二维主成分空间，结果如下图：

image comes from https://projector.tensorflow.org/

从图可以看出：

不同的数字还是有聚集效果，比如1集中在左边，0集中在右边;
不标记颜色和数字的话，所有数字集中在一起。

使用t-SNE(一种非线性降维方法)基本能把这些数字分开^t-SNE。

^智能时代. 《智能时代》, 吴军. 第三章: 思维的革命, 小节: 熵,一种新的世界观。 ↩

^PRML. Pattern Recognition and Machine Learning, Bishop. ↩

^t-SNE. Visualizing High-Dimensional Data Using t-SNE, 2008. First author’s blog: https://lvdmaaten.github.io/tsne/ ↩