Machine Learning for Graphs and Sequential Data

Sequential Data / Temporal Data

针对Sequential Data我们一般会这样划分机器学习的数据集，只用未来的数据进行测试：

Autoregressive Models

Motivation & Definitions

Motivation: 很多时间序列都有这种特点：

今天气温和昨天气温接近
今天销量会受前几天销量影响
今天利率、通胀、库存、流量，通常和前几天的数据有关系

所以就有了以下定义的 autoregressive model：

Definition: An autoregressive model AR(p) of order $p$ is defined as

$X_t = c + \sum_{i=1}^{p} \varphi_i X_{t-i} + \varepsilon_t$

where $\varphi_1, \ldots, \varphi_p$ are the parameters, $c$ is a constant, and $\varepsilon_t \sim N(0,\sigma)$ is white noise. The variable $X_{t-i}$ is the lagged value of $X$ at time $t-i$.

也可以写成：

$P(X_t \mid X_{t-1}, \ldots, X_{t-p}) \sim N\left(c + \sum_{i=1}^{p} \varphi_i X_{t-i}, \sigma\right)$

注意，随着t增加，$\varphi_i$不变，但是会跟着一起”平移“：

$X_t = c + \varphi_1 X_{t-1} + \varphi_2 X_{t-2} + \cdots + \varphi_p X_{t-p} + \varepsilon_t$ $X_{t+1} = c + \varphi_1 X_t + \varphi_2 X_{t-1} + \cdots + \varphi_p X_{t+1-p} + \varepsilon_{t+1}$

对应的Graphical Model：

Definition:

The mean function of an AR model is $\mu(t) = E[X_t]$.
By default, it depends on $t$.
The autocovariance is $\gamma(t,i) = \operatorname{Cov}(X_t, X_{t-i})$.
By default, it depends on $t$ and $i$.
The autocovariance function can be normalized to give the Pearson autocorrelation function
$\rho(t,i) = \frac{\operatorname{Cov}(X_t, X_{t-i})}{\sqrt{\operatorname{Var}(X_t)}\sqrt{\operatorname{Var}(X_{t-i})}}.$
It lies in $[-1,1]$.

参数学习（Parameter learning）

先假设$c=0$，不难发现对 $t = p, p+1, p+2, \ldots, T$，分别有：：

$\begin{align*} X_p &= \varphi_1 X_{p-1} + \varphi_2 X_{p-2} + \cdots + \varphi_p X_0 + \varepsilon_p \\ X_{p+1} &= \varphi_1 X_p + \varphi_2 X_{p-1} + \cdots + \varphi_p X_1 + \varepsilon_{p+1} \\ X_{p+2} &= \varphi_1 X_{p+1} + \varphi_2 X_p + \cdots + \varphi_p X_2 + \varepsilon_{p+2} \\ &\vdots \end{align*}$

写成矩阵的形式就是：

$\begin{bmatrix} X_p \\ X_{p+1} \\ X_{p+2} \\ \vdots \\ X_T \end{bmatrix} = \begin{bmatrix} X_{p-1} & X_{p-2} & \cdots & X_0 \\ X_p & X_{p-1} & \cdots & X_1 \\ X_{p+1} & X_p & \cdots & X_2 \\ \vdots & \vdots & \ddots & \vdots \\ X_{T-1} & X_{T-2} & \cdots & X_{T-p} \end{bmatrix} \begin{bmatrix} \varphi_1 \\ \varphi_2 \\ \vdots \\ \varphi_p \end{bmatrix} + \begin{bmatrix} \varepsilon_p \\ \varepsilon_{p+1} \\ \varepsilon_{p+2} \\ \vdots \\ \varepsilon_T \end{bmatrix}$

记成：

$y = X\varphi + \varepsilon$

用OLS（ordinary least squares）的公式可得：

$\hat{\varphi} = (X^TX)^{-1}X^Ty$

如果$c$不等于0的话，把$X$和$\varphi$调整成：

$\begin{bmatrix} 1 & X_{p-1} & X_{p-2} & \cdots & X_0 \\ 1 & X_p & X_{p-1} & \cdots & X_1 \\ 1 & X_{p+1} & X_p & \cdots & X_2 \\ \vdots & \vdots & \ddots & \vdots \\ 1 & X_{T-1} & X_{T-2} & \cdots & X_{T-p} \end{bmatrix} \begin{bmatrix} c\\ \varphi_1 \\ \varphi_2 \\ \vdots \\ \varphi_p \end{bmatrix}$

再套用公式即可。

Stationary

Definition:

A process is said to be stationary if

1.$E[X_t] = E[X_{t-i}] = \mu,\ \forall\ t,\forall\ i$（也就是$E[X_t] \equiv \mu \ \forall t$）

2. $\text{Cov}(X_t, X_{t-i}) = \gamma_i,\ \forall\ t,\forall\ i$

3. $E\left[|X_t|^2\right] < \infty,\ \forall\ t$

Theorem: Consider the stationary AR$(p)$ process

$X_t=c+\sum_{i=1}^{p}\phi_i X_{t-i}+\epsilon_t$

with $\epsilon \sim \mathcal{N}(0,\sigma^2)$. Then:

$\begin{align*}\mu&=\frac{c}{1-\sum_{i=1}^{p}\phi_i},\quad \text{for all } t \\\gamma_0&=\sum_{j=1}^{p}\phi_j\gamma_{-j}+\sigma^2 \\\gamma_i&=\sum_{j=1}^{p}\phi_j\gamma_{i-j},\quad \text{for all } t,\ i\in[1,p]\end{align*}$

Proof：

$\begin{align*} \mathbb{E}[X_t] &= c+\sum_{i=1}^{p}\varphi_i\mathbb{E}[X_{t-i}] +\mathbb{E}[\epsilon_t] \\ \Rightarrow \mu &= c+\sum_{i=1}^{p}\varphi_i\mu+0 \\ \Rightarrow \mu &= \frac{c}{1-\sum_{i=1}^{p}\varphi_i} \end{align*}$

2,3.

$\begin{align*} \operatorname{Cov}(X_t,X_{t-i}) &= \operatorname{Cov}(c,X_{t-i}) + \sum_{j=1}^{p}\phi_j \operatorname{Cov}(X_{t-j},X_{t-i}) + \operatorname{Cov}(\epsilon_t,X_{t-i}) \\ \Rightarrow \gamma_i &= \sum_{j=1}^{p}\phi_j\gamma_{i-j} + \operatorname{Cov}(\epsilon_t,X_{t-i}) \end{align*}$

因为$\epsilon_t,\epsilon_{t-i}$之间是independent，所以$\epsilon_t$也和之前的$X_{t-i}$是independent的：

$\operatorname{Cov}(\epsilon_t,X_{t-i}) = \begin{cases} \sigma^2, & \text{if } i=0, \\ 0, & \text{otherwise}. \end{cases}$

从而得到

$\begin{align*} \gamma_0 &= \sum_{j=1}^{p}\phi_j\gamma_{-j} + \sigma^2 \\ \gamma_i &= \sum_{j=1}^{p}\phi_j\gamma_{i-j}, \quad \text{for all } t,\ i\in[1,p] \end{align*}$

$\square$

在满足stationary的条件下，我们可以用Yule-Walker equations来学习参数：

Yule-Walker equations：

$\begin{bmatrix} \gamma_1 \\ \gamma_2 \\ \vdots \\ \gamma_{p-1} \\ \gamma_p \end{bmatrix} = \begin{bmatrix} \gamma_0 & \gamma_{-1} & \cdots & \gamma_{2-p} & \gamma_{1-p} \\ \gamma_1 & \gamma_0 & \cdots & \gamma_{3-p} & \gamma_{2-p} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ \gamma_{p-2} & \gamma_{p-3} & \cdots & \gamma_0 & \gamma_{-1} \\ \gamma_{p-1} & \gamma_{p-2} & \cdots & \gamma_1 & \gamma_0 \end{bmatrix} \begin{bmatrix} \varphi_1 \\ \varphi_2 \\ \vdots \\ \varphi_{p-1} \\ \varphi_p \end{bmatrix}$

因为

$\gamma_{-i} = \text{Cov}(X_t,X_{t+i}) = \text{Cov}(X_{t+i},X_t) = \gamma_{i}$

所以这里的Yule-Walker矩阵也可以写成对称的形式：

$\begin{bmatrix} \gamma_0 & \gamma_{1} & \cdots & \gamma_{p-2} & \gamma_{p-1} \\ \gamma_1 & \gamma_0 & \cdots & \gamma_{p-3} & \gamma_{p-2} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ \gamma_{p-2} & \gamma_{p-3} & \cdots & \gamma_0 & \gamma_{1} \\ \gamma_{p-1} & \gamma_{p-2} & \cdots & \gamma_1 & \gamma_0 \end{bmatrix}$

具体流程：

1. 计算协方差$\gamma_0,…,\gamma_p$：

先估计均值：

$\bar{x} = \frac{1}{T} \sum_{t=1}^{T} x_t$

然后用样本自协方差估计 moments：

$\hat{\gamma}_k = \frac{1}{T} \sum_{t=k+1}^{T} (x_t-\bar{x})(x_{t-k}-\bar{x}), \quad k=0,1,\ldots,p$

2. 计算Yule_Walker矩阵的逆，解得$\varphi_1,…,\varphi_p$ （解方程）

3. 利用$\gamma_0$的公式来estimate $\sigma$：

$\gamma_0 = \sum_{j=1}^{p}\phi_j\gamma_{-j} + \sigma^2$

注意到其实Yule-Walker矩阵本质上是一个协方差矩阵，假设

$Z_t = \begin{bmatrix} X_t \\ X_{t-1} \\ \vdots \\ X_{t-p+1} \end{bmatrix}$

那么

$\text{Cov}(Z_t) = \begin{bmatrix} \gamma_0 & \gamma_{1} & \cdots & \gamma_{p-1} \\ \gamma_1 & \gamma_0 & \cdots & \gamma_{p-2} \\ \vdots & \vdots & \ddots & \vdots \\ \gamma_{p-1} & \gamma_{p-2} & \cdots & \gamma_0 \end{bmatrix}$

Markov Chains

我们假设r.v $X_t$和time index都是离散的。

Definition: A Markov Chain is a sequence of r.v. $X_1, X_2, \ldots, X_T$ which fulfills the Markov property:

$P(X_t \mid X_1, \ldots, X_{t-1}) = P(X_t \mid X_{t-1})$

也就是说，当前状态出现的概率只取决于上一个状态。

Joint distribution则是：

$P(X_1 = i_1, \ldots, X_T = i_T) = P(X_1 = i_1)\prod_{t=1}^{T-1} P(X_{t+1} = i_{t+1} \mid X_t = i_t)$

General的情况下，每个r.v的分布可能是不同的，也就是说

$P(X_1 = i) = \pi_i \quad \text{and} \quad P(X_{t+1} = j \mid X_t = i) = A_{ij}^{(t+1)}$

here $\pi \in \mathbb{R}^K$ is a prior probability on the initial state, and $A^{(t)} \in \mathbb{R}^{K \times K}$ are the transition matrices.

但是我们这里考虑简单点的情况：time-homogeneous / stationary Markov Chain

$P(X_1 = i) = \pi_i \quad \text{and} \quad P(X_{t+1} = j \mid X_t = i) = A_{ij}$

也就是$A^{(t)} \equiv A \ \forall t$。

Theorem:

$P(X_{t+n}=j \mid X_t=i) = A^n_{ij}$

Theorem:

$P(X_t=j) = (\pi A^{(t-1)})_j$

参数学习（Parameter learning）

核心思路就是用MLE（maximum likelihood estimation）。

假设有$N$条观测序列：

$\mathcal{D} = \{ x^{(1)}_{1:T_1}, x^{(2)}_{1:T_2}, ... , x^{(N)}_{1:T_N} \}$

每条的长度是$T_i$。比如说

$\begin{align*} x^{(1)}_{1:3} &= [1,3,2] \\ x^{(2)}_{1:1} &= [3] \\ x^{(3)}_{1:4} &= [1,1,3,2] \\ \end{align*}$

对于每条观测序列都有：

$\begin{align*} P(x_{1:T}) &= P(x_1)\prod_{t=1}^{T-1} P(x_{t+1}\mid x_t) \\ & = \pi_{x_1}\prod_{t=1}^{T-1} A_{x_t x_{t+1}} \end{align*}$

所以

$\begin{align*} P(\mathcal{D}) &= \prod_{n=1}^{N} P\!\left(x_{1:T_n}^{(n)}\right) \\ & = \prod_{n=1}^{N} \left( \pi_{x_1^{(n)}} \prod_{t=1}^{T_n-1} A_{x_t^{(n)},x_{t+1}^{(n)}} \right) \end{align*}$

统计初始状态次数（多少条序列是从$i$开始的）以及状态转移次数（$i$到$j$的次数）：

$c(i) = \# (X_1 = i), \ c(i,j) = \# (X_t = i,X_{t+1} = j)$

于是

$\log P(\mathcal{D}) = \sum_i c(i)\log \pi_i + \sum_i \sum_j c(i,j)\log A_{ij}$

在满足$\sum_k \pi_k = 1$和$\sum_j A_{ij} = 1$的情况下maximize $\log P(\mathcal{D})$就会得到：

$\hat{\pi}_i = \frac{c(i)}{\sum_{i'} c(i')}$ $\hat{A}_{ij} = \frac{c(i,j)}{\sum_{j'} c(i,j')}$

所以$\hat{A}$适合一行一行算。

例子：

$\begin{align*} X^{(1)} &= [1,3,2] \\ X^{(2)} &= [3] \\ X^{(3)} &= [1,1,3,2] \\ \end{align*}$

统计一下：

$c(1)=2,c(2)=0,c(3)=1$ $c(1,1)=1, c(1,3) = 2, c(3,2)=2$

所以可以得到

$\pi = \begin{pmatrix} \frac{2}{3} \\ 0 \\ \frac{1}{3} \end{pmatrix}$ $A = \begin{pmatrix} \frac{1}{3} & 0 & \frac{2}{3}\\ 0 & 0 & 0\\ 0 & 1 & 0 \end{pmatrix}$

Hidden Markov Models

现实中很多系统的真实状态并不能直接被看到，我们只能通过观测来间接感知。

所以我们需要考虑当序列数据的真实状态看不见，只能看到带噪声的观测时，怎么建模、推断和学习参数。

比如说我们考虑一个观测飞机位置的例子：

真正的状态是位置、速度等物理量，是隐藏变量
传感器读到的是有噪声的位置，是观测变量

真实状态序列可能满足 Markov 性，但观测序列本身未必满足，所以需要 HMM。

Definition: A Hidden Markov Model (HMM) is composed of a sequence of hidden/latent variables $[Z_1,\ldots,Z_T]$ and a sequence of observed variables $[X_1,\ldots,X_T]$ such that:

The r.v. $Z_1,\ldots,Z_T$ satisfy the Markov property:

$P(Z_{t+1}\mid Z_t,Z_{t-1},\ldots,Z_1)=P(Z_{t+1}\mid Z_t)$

where $P(Z_{t+1}\mid Z_t)$ are the transition probabilities.

Distribution of $X_t$ depends only on $Z_t$:

$P(X_{t+1}\mid Z_1,\ldots,Z_T,X_1,\ldots,X_T)=P(X_{t+1}\mid Z_{t+1})$

where $P(X_{t+1}\mid Z_{t+1})$ are the emission probabilities.

也就是说当前隐藏状态只依赖上一个隐藏状态，并且当前观测状态只取决于当前的隐藏状态。

和之前一样，可以这样建模：

$\begin{align*} P(Z_1=i)&=\pi_i \\ P(Z_{t+1}=j \mid Z_t=i)&=A_{ij}^{(t+1)} \\ P(X_{t+1}=j \mid Z_{t+1}=i)&=B_{ij}^{(t+1)} \\ \end{align*}$

但为了降低参数数量，我们还是假设所有时间点共享一样的分布$A,B$：

在HMM中，我们一般有2种任务（都是假设已知观测数据$X_{1:T}$）：

Inference：参数（$(\pi, A, B)$）已知，推导隐藏状态
Parameter Learning：隐藏状态未知，学习模型参数

Inference

而Inference一般也分3种：

Filtering: computes $P(Z_t \mid X_{1:t} )$ online

只拥有到目前为止的所有观测信息，也就是及时判断。
Smoothing: computes $P(Z_t \mid X_{1:T} )$ offline

拥有全部的观测信息，也就是事后回顾。
MAP inference: computes $\text{arg } \underset{Z_{1:T}}{\text{max}}P(Z_{1:T} \mid X_{1:T} )$

Filtering

computes $P(Z_t \mid X_{1:t} )$.

Forward Algorithm:

记

$\alpha_t(k) := P(Z_t = k, X_{1:t}) \quad \text{and} \quad \alpha_t := \begin{bmatrix} \alpha_t(1)\\ \vdots \\ \alpha_t(K)\\ \end{bmatrix}$

根据Bayes rule：

$\begin{align*} P(Z_t=k \mid X_{1:t}) &= \frac{P(Z_t=k, X_{1:t})} {\sum_{j=1}^{K} P(Z_t=j, X_{1:t})} \\ &= \frac{\alpha_t(k)} {\mathrm{sum}(\boldsymbol{\alpha}_t)} \end{align*}$

Forward Algorithm主要就分2步：

计算$\alpha_1$（initialisation）：
$\begin{align*} \alpha_1(k) &= P(Z_1=k, X_1) \\ &= P(Z_1=k)P(X_1 \mid Z_1=k) \\ &= \pi_k B_{k x_1} \end{align*}$
给定$\alpha_t$，计算$\alpha_t$（recursion）：
$\begin{align*} \alpha_{t+1}(k) &= P(Z_{t+1}=k, X_{1:t+1}) \\ &= P(X_{t+1}\mid Z_{t+1}=k, X_{1:t}) P(Z_{t+1}=k, X_{1:t}) \\ &= P(X_{t+1}\mid Z_{t+1}=k) \sum_{j=1}^{K} P(Z_{t+1}=k, Z_t=j, X_{1:t}) \\ &= P(X_{t+1}\mid Z_{t+1}=k) \sum_{j=1}^{K} P(Z_{t+1}=k \mid Z_t=j, X_{1:t}) P(Z_t=j, X_{1:t}) \\ &= B_{k(x_{t+1})} \sum_{j=1}^{K} A_{jk}\alpha_t(j) \end{align*}$
写成矩阵的话就是：
$\boldsymbol{\alpha}_{t+1} = \mathbf{B}_{:(x_{t+1})} \odot \left(A^T\boldsymbol{\alpha}_t\right)$ $\begin{bmatrix} \alpha_{t+1}(1) \\ \alpha_{t+1}(2) \\ \vdots \\ \alpha_{t+1}(K) \end{bmatrix} = \begin{bmatrix} B(1,x_{t+1}) \\ B(2,x_{t+1}) \\ \vdots \\ B(K,x_{t+1}) \end{bmatrix} \odot \left( A^T \begin{bmatrix} \alpha_t(1) \\ \alpha_t(2) \\ \vdots \\ \alpha_t(K) \end{bmatrix} \right)$
（其中$\odot$ denotes the Hadamard product：$(A \odot B)_{ij} = (A)_{ij}(B)_{ij}.$）

计算$\alpha_{1:T}$的复杂度为$O(TK^2)$。

Smoothing

computes $P(Z_t \mid X_{1:T} )$ .

Forward-Backward Algorithm：

在之前forward的基础上再定义一个$\beta$：

$\beta_t(k) := P(X_{t+1:T}\mid Z_t=k) \quad \text{and} \quad \boldsymbol{\beta}_t = \begin{bmatrix} \beta_t(1) \\ \vdots \\ \beta_t(K) \end{bmatrix}$

那么根据Bayes就有：

$\begin{align*} P(Z_t=k \mid X_{1:T}) &= \frac{ P(Z_t=k, X_{1:t})P(X_{t+1:T}\mid Z_t=k) }{ \sum_{j=1}^{K} P(Z_t=j, X_{1:T}) } \\ &\propto \alpha_t(k)\beta_t(k) \end{align*}$

Backward Algorithm也分2步：

计算$\beta_T$（initialisation）：
$\beta_T(k)= 1$
因为在$T$的时候：
$P(Z_t=T \mid X_{1:T})\propto \alpha_T(k)\beta_T(k)$
$\alpha_T(k)$已经捕获了所有已知信息了，所以$\beta_T(k)$就需要是个常量。
给定$\beta_{t+1}$，计算$\beta_t$（recursion）：
$\begin{align*} \beta_t(j) &= P(X_{t+1:T}\mid Z_t=j) \\ &= \sum_{k=1}^{K} P(X_{t+1:T}, Z_{t+1}=k \mid Z_t=j) \\ &= \sum_{k=1}^{K} P(X_{t+1}, Z_{t+1}=k \mid Z_t=j) P(X_{t+2:T}\mid Z_t=j, X_{t+1}, Z_{t+1}=k) \\ &= \sum_{k=1}^{K} P(Z_{t+1}=k\mid Z_t=j) P(X_{t+1}\mid Z_{t+1}=k, Z_t=j) P(X_{t+2:T}\mid Z_{t+1}=k) \\ &= \sum_{k=1}^{K} A_{jk}B_{k x_{t+1}}\beta_{t+1}(k) \end{align*}$
写成矩阵的话就是：
$\boldsymbol{\beta}_t = A\left( \mathbf{B}_{:(x_{t+1})} \odot \boldsymbol{\beta}_{t+1} \right)$ $\begin{bmatrix} \beta_t(1) \\ \beta_t(2) \\ \vdots \\ \beta_t(K) \end{bmatrix} = A \left( \begin{bmatrix} B(1,x_{t+1}) \\ B(2,x_{t+1}) \\ \vdots \\ B(K,x_{t+1}) \end{bmatrix} \odot \begin{bmatrix} \beta_{t+1}(1) \\ \beta_{t+1}(2) \\ \vdots \\ \beta_{t+1}(K) \end{bmatrix} \right)$

计算$\beta_{1:T}$的复杂度也为$O(TK^2)$。

最后计算

$\gamma_t(k) := P(Z_t=k \mid X_{1:T}) = \frac{\alpha_t(k)\beta_t(k)} {\sum_s \alpha_t(s)\beta_t(s)}$

即可。

MAP inference

Viterbi algorithm：

注意到：

$\begin{align*} \arg\max_Z P(Z_{1:T}\mid X_{1:T}) &= \arg\max_Z \log\left[P(Z_{1:T},X_{1:T})\right] \\ &= \arg\max_Z \left(\log\left[P(Z_1)P(X_1\mid Z_1)\right] \ + \sum_{t=2}^{T} \log\left[ P(Z_t\mid Z_{t-1})P(X_t\mid Z_t) \right] \right) \end{align*}$

注意到每一项都只取决于$Z_t$和$Z_{t-1}$，所以我们可以把它转成一个图论里的最短路径问题：

其中：

和start相连的边的权重：$-\log\left[ \Pr(Z_1=j)\Pr(X_1\mid Z_1=j) \right]$
中间层：$-\log\left[\Pr(Z_t=j\mid Z_{t-1}=i)\Pr(X_t\mid Z_t=j)\right]$
和start相连的边：$0$

复杂度依旧是$O(TK^2)$。

Parameter Learning

$X_{1:T}$已知，$Z_{1:T}$未知。

目标：求

$\begin{align*} \underset{\theta}{\text{max}} \log P(X\mid \theta) = \underset{\theta}{\text{max}} \log \sum_Z P(X,Z \mid \theta) \end{align*}$

也就是说虽然 $X$ 观测到了,但背后的状态序列 $Z_{1:T}$ 没看到,所以必须把所有可能的隐藏状态序列都加起来。但隐藏状态序列的可能数是指数级的，所以不能直接做。

但注意到这和GMM问题很像，所以我们考虑用EM算法求近似解。

Baum-Welch algorithm：

E step: evaluate the posterior：
$P(Z\mid X,\theta^{old}) = P(Z_{1:T}\mid X_{1:T},\theta^{old}) = \frac{ P(X_{1:T}\mid Z_{1:T},\theta^{old}) P(Z_{1:T}\mid \theta^{old}) }{ \sum_{Z_{1:T}} P(X_{1:T}\mid Z_{1:T},\theta^{old}) P(Z_{1:T}\mid \theta^{old}) }$
其中
$P(X_{1:T}\mid Z_{1:T},\theta^{old}) P(Z_{1:T}\mid \theta^{old}) = \prod_{t=1}^{T} P(X_t\mid Z_t,\theta^{old}) \cdot \prod_{t=2}^{T} P(Z_t\mid Z_{t-1},\theta^{old}) \cdot P(Z_1\mid \theta^{old})$
M step: maximize the expected joint log-likelihood
$\theta^{new} = \arg\max_{\theta} \mathbb{E}_{P(Z\mid X,\theta^{old})} \left[ \log p(X,Z\mid \theta) \right]$
其中
$\begin{align*} \mathbb{E}_{P(Z_{1:T}\mid X_{1:T},\theta^{old})} \left[ \log P(X_{1:T},Z_{1:T}\mid \theta) \right] &= \sum_k P(Z_1=k\mid X_{1:T},\theta^{old}) \log(\pi_k) \\ &\quad+ \sum_{i,j}\sum_t P(Z_t=i,Z_{t+1}=j\mid X_{1:T},\theta^{old}) \log(A_{ij}) \\ &\quad+ \sum_i\sum_t P(Z_t=i\mid X_{1:T},\theta^{old}) \mathbb{I}(x_t=j) \log(B_{ij}) \end{align*}$

Continuous Data

我们之前假设的是离散的情况：discrete time $t \in \{1,2,\ldots,T\}$ and discrete r.v.$Z_t \in \{1,2,\ldots,K\}$, $X_t \in \{1,2,\ldots,K’\}$

$\begin{align*} P(Z_1=i)&=\pi_i \\ P(Z_{t+1}=j\mid Z_t=i)&=A_{ij} \\ P(X_{t+1}=j\mid Z_{t+1}=i)&=B_{ij} \\ \end{align*}$

我们现在关注discrete time $t \in \{1,2,\ldots,T\}$, discrete r.v. $Z_t \in \{1,2,\ldots,K\}$, and continuous $X_t \in \mathbb{R}^d$:

$\begin{align*} P(Z_1=i)&=\pi_i \\ P(Z_{t+1}=j\mid Z_t=i)&=A_{ij} \\ P(X_{t+1}=x\mid Z_{t+1}=i) &= N(x\mid \mu_i,\sigma_i^2 \cdot I) \\ \end{align*}$

Inference部分基本不变，只需要把所有计算

$P(X_{t+1}=j\mid Z_{t+1}=i)=B_{ij}$

的部分改成

$P(X_{t+1}=x\mid Z_{t+1}=i) = N(x\mid \mu_i,\sigma_i^2 \cdot I)$

即可。

而Parameter Learning的部分，E-Step没有变化，M-step需要一点调整：

$\begin{align*} \mathbb{E}_{P(Z_{1:T}\mid X_{1:T},\theta^{old})} \left[ \log P(X_{1:T},Z_{1:T}\mid \theta) \right] &= \sum_k P(Z_1=k\mid X_{1:T},\theta^{old}) \log(\pi_k) \\ &\quad+ \sum_{i,j}\sum_t P(Z_t=i,Z_{t+1}=j\mid X_{1:T},\theta^{old}) \log(A_{ij}) \\ &\quad+ \sum_i\sum_t P(Z_t=i\mid X_{1:T},\theta^{old}) \log\left( N(X_t\mid \mu_i,\sigma_i^2 \cdot I) \right) \end{align*}$

而它的参数更新的公式：

$\mu_i^{new} = \frac{ \sum_{t=1}^{T}\gamma_t(i)X_t }{ \sum_{t=1}^{T}\gamma_t(i) }$ $\sigma_i^{new} = \frac{ \sum_{t=1}^{T} \gamma_t(i) (\mu_i^{new}-x_t)(\mu_i^{new}-x_t)^T }{ \sum_{t=1}^{T}\gamma_t(i) } \quad \text{where } \gamma_t(i) := P(Z_t=i\mid X_{1:T},\theta^{old})$

和GMM里更新高斯参数的形式是等价/一样的。

唯一的区别就是：

在 GMM 里，不同观测样本$X_t$是独立的（independent）
在 HMM 里，观测 $X_t$ 在给定隐藏状态之后是条件独立的（conditional independent given $Z$）

Neural Network Approches

Word Vectors

由于神经网络不能直接处理单词，所以我们需要先把词表示成向量，方便后续计算处理。

最简单的想法就是one-hot encoding：假设有个单词数量为n的单词表，把每个单词转成一个维度为n的向量（其中有且仅有一项等于1）。但它有三个问题：维度高、稀疏、假设词之间相互独立。

于是就引出了 word vectors / embeddings：相似词的向量应该接近，比如 “hotel” 和 “motel”，”man”和”woman”，”walking”和”walked”。这基于分布式假设：出现在相似上下文中的词，通常有相似含义。

有2种常见的获取word vectors的方法：

1. Co-occurence Matrix + SVD：

这个方法偏统计。它的初始想法很直观：一个词的含义可以通过“它经常和哪些词一起出现”来表示。

先统计一个词附近出现哪些词，得到一个co-occurrence matrix。然后再用SVD把高维稀疏向量压缩成低维稠密向量。

优点是相似的单词会得到相似的向量。

缺点是计算慢，而且新词不好加入。

2. Word2Vec

这种方法偏预测。

先设计一个预测任务，然后用神经网络学习词向量。

来看2个经典模型：

1. CBOW (Continuous Bag-of-Words)：

给定一个词周围的上下文词（window），预测中间词。

但是这种方法对稀有词不太好，因为稀有词出现次数非常少，模型很少有机会把它作为要预测的目标词，所以不容易学好它的向量。

2. Skip-gram：

给定中心词，预测它周围的上下文词。

它对稀有词更友好（可以这样理解：如果一个 rare word 出现了，模型会强迫自己用这个 rare word 去预测它的上下文。这样模型必须理解这个 rare word 通常出现在什么环境中，因此更容易学到它的含义。），但训练起来也会更慢。

来看一下具体的结构：

假设词表大小是$N$，并且每个词一开始用one-hot vector来表示。

Embedding则是通过一个矩阵$U \in \mathbb{R}^{N \times D}$来实现的。由于所有输入都是one-hot vector，所有相当于$U$的每一行都是对应一个单词的embedding结果的。

在此基础上还需要一个矩阵$V \in \mathbb{R}^{D \times N}$，用来把中心词 embedding 转换成对整个词表的预测分数。最后再套上softmax就会得到一个上下文单词的概率分布。

也就是整体流程是这样的：

$\text{softmax}(x_i^T \cdot U \cdot V) = \text{softmax}(u_i V)$

假设$S=[x_{i-l},\ldots,x_{i-1},x_{i+1},\ldots,x_{i+l}]$是当前$x_i$的大小为$l$的window，$\theta = (U,V)$为模型参数，那么我们的训练目标就是最大化这个期望值：

$\max_{\theta}\mathbb{E}\left[P(S\mid x_i,\theta)\right] = \min_{\theta}\left(-\mathbb{E}\left[P(S\mid x_i,\theta)\right]\right)$

其中

$\begin{align*} P(S\mid x_i,\theta) = \prod_{x_k\in S} P(x_k\mid x_i,\theta) = \prod_{x_k\in S} \operatorname{softmax}(u_i V)_k \end{align*}$

但是注意到，softmax的形式是：

$P(x_k \mid x_i,\theta) = \frac{ \exp(u_i v_k^T) }{ \sum_{l=1}^{N}\exp(u_i v_l^T) }$

为了计算这个分母我们需要遍历整个单词表，当单词表很大的时候会非常低效。

所以可以用negative sampling技巧：我们不再把问题看成从整个词表中预测正确上下文词，而是把它变成一个二分类问题：判断一个词对是不是真实的上下文关系。

利用（最小化）这个Loss来训练模型：

$L = \log\left(P(x_p\mid x_i,\theta)\right) + \log\left(1-P(x_n\mid x_i,\theta)\right)$

并且用sigmoid代替softmax：

$P(x_k\mid x_i,\theta) = \operatorname{sigmoid}(u_i v_k^T)$

如果两个词经常一起出现，它们的向量点积应该大，sigmoid 后接近 1。
如果两个词很少一起出现，它们的向量点积应该小，sigmoid 后接近 0。

RNNs

我们刚才看了如何处理单个单词。但我们该怎么处理一整句话、一个音频片段、一个时间序列呢？尤其是它们的长度都各不相同。

比如说当我们希望它接收输入：

$x=[I,like,dogs]$

，并且输出

$y=[PRON,VERB,NOUN]$

答案就是Recurrent Neural Networks (RNNs)。它的基本想法是用一个 hidden state 来总结到目前为止读过的序列信息。

RNN可以处理多种任务：

来看正式的定义：

假设输入是$\{x^{(1)},\ldots,x^{(T)}\}$，并且输出$\{y^{(1)},\ldots,y^{(T)}\}$。我们现在希望知道

$P(y^{(t)}\mid x^{(1)},\ldots,x^{(t)})$

的概率分布。

所以我们用 $h^{(t-1)}$表示/记录$\{x^{(1)},\ldots,x^{(t-1)}\}$的信息。

具体的更新公式：

$\begin{align*} z^{(t)} &= W h^{(t-1)} + Ux^{(t)} + b \\ h^{(t)} &= \tanh(z^{(t)}) \\ o^{(t)} &= Vh^{(t)} \\ \hat{y}^{(t)} &= \operatorname{softmax}(o^{(t)}) \end{align*}$

训练的时候关注negative log-likelihood：

$\begin{align*} L &= -\log \prod_t p_{\mathrm{model}} \left( y^{(t)}\mid x^{(1)},\ldots,x^{(t)} \right) \\ &= -\sum_t \log p_{\mathrm{model}} \left( y^{(t)}\mid x^{(1)},\ldots,x^{(t)} \right) \\ &= -\sum_t L^{(t)} \end{align*}$

当我们将RNN展开（unrolling）：

就会发现它像一个很深的神经网络，所以可以用反向传播训练。

假设真实标签$y^{(t)}$用的是one-hot label（即其中只有一位等于1，其余位均等于0），$\hat{y}_c^{(t)}$为每轮训练时得到的softmax之后的分布，那么$L^{(t)}$就可以写成

$\begin{align*} L^{(t)} = \text{log} \prod_{c=1}^{C} \left(\hat{y}_c^{(t)}\right)^{y_c^{(t)}} = \sum_{c=1}^{C} y_c^{(t)} \log \hat{y}_c^{(t)} \end{align*}$

来计算它的导数。

先看$\frac{\partial L^{(t)}}{\partial o_i^{(t)}}$：

$\begin{align*} L^{(t)} &= \sum_c y_c^{(t)} \log\left( \frac{ e^{o_c^{(t)}}}{ \sum_i e^{o_i^{(t)}} } \right) \\ &= \sum_c y_c^{(t)} \left( o_c^{(t)} - \log\left(\sum_i e^{o_i^{(t)}}\right) \right) \\ &= \sum_c y_c^{(t)}o_c^{(t)} - \sum_c y_c^{(t)} \log\left(\sum_i e^{o_i^{(t)}}\right) \\ &= \sum_c y_c^{(t)}o_c^{(t)} - \log\left(\sum_i e^{o_i^{(t)}}\right) \end{align*}$

所以

$\frac{\partial L^{(t)}}{\partial o_i^{(t)}} = y_i^{(t)} -\hat{y}_i^{(t)}$ $\frac{\partial L^{(t)}}{\partial o^{(t)}} = y^{(t)} -\hat{y}^{(t)}$

再来看$\frac{\partial L^{(t)}}{\partial V^{(t)}}$

$\begin{align*} \frac{\partial L^{(t)}}{\partial V^{(t)}} &= \frac{\partial L^{(t)}}{\partial o^{(t)}} \cdot\frac{\partial o^{(t)}}{\partial V^{(t)}} \\ &= (y^{(t)} -\hat{y}^{(t)}) \cdot (h^{(t)})^T \end{align*}$ $\begin{align*} \frac{\partial L}{\partial V} &= - \sum_t (y^{(t)} -\hat{y}^{(t)}) \cdot (h^{(t)})^T \\ &= \sum_t (\hat{y}^{(t)} - y^{(t)}) \cdot (h^{(t)})^T \end{align*}$

类似的有：

$\begin{align*} \frac{\partial L^{(t)}}{\partial W^{(t)}} &= \frac{\partial L^{(t)}}{\partial h^{(t)}} \cdot \frac{\partial h^{(t)}}{\partial z^{(t)}} \cdot \frac{\partial z^{(t)}}{\partial W^{(t)}}\\ &= \frac{\partial L}{\partial h^{(t)}} \text{diag}\left(1-\left(h^{(t)}\right)^2\right) \left(h^{(t-1)}\right)^T \end{align*}$ $\frac{\partial L}{\partial W} = \sum_t \frac{\partial L}{\partial h^{(t)}} \text{diag}\left(1-\left(h^{(t)}\right)^2\right) \left(h^{(t-1)}\right)^T$ $\frac{\partial L}{\partial U} = \sum_t \frac{\partial L}{\partial h^{(t)}} \text{diag}\left(1-\left(h^{(t)}\right)^2\right) \left(x^{(t)}\right)^T$

但注意到里面用到的$h^{(t)}$是递归定义的：

$h^{(t)} = \tanh\left(Wh^{(t-1)} + Ux_t + b\right)$

所以当计算$\frac{\partial L}{\partial h^{(t)}}$时：

$\begin{align*} \frac{\partial L}{\partial h^{(t)}} &= \frac{\partial L^{(t)}}{\partial h^{(t)}} + \frac{\partial L^{(t+1)}}{\partial h^{(t)}} + \cdots + \frac{\partial L^{(T)}}{\partial h^{(t)}} \\ &= \frac{\partial L^{(t)}}{\partial h^{(t)}} + \frac{\partial L^{(t+1)}}{\partial h^{(t+1)}} \cdot \frac{\partial h^{(t+1)}}{\partial h^{(t)}} + \cdots + \frac{\partial L^{(T)}}{\partial h^{(T)}} \cdot \frac{\partial h^{(T)}}{\partial h^{(t)}} \end{align*}$

其中

$\frac{\partial L^{(i)}}{\partial h^{(i)}} = \left(\hat{y}^{(i)} - y^{(i)}\right)^T V$ $\frac{\partial h^{(e)}}{\partial h^{(s)}} = \prod_{i=s}^{e} \frac{\partial h^{(i+1)}}{\partial h^{(i)}} = \prod_{i=s}^{e} W\left(1-\left(h^{(i)}\right)^2\right)$

$W$的连乘会导致这个梯度消失或者爆炸（vanish/explode），进而会导致RNN无法有效地更新/学习到信息。

为了解决这个问题我们需要完善目前的架构。

GRU (Gated Recurrent Unit)

核心思想：不是每一步都完全覆盖 hidden state，而是引入gate来控制保留多少旧信息以及接受多少新信息。

（有点类似residual NN的idea。）

更新公式：

$\begin{align*} z^{(t)}&= \sigma\left(W_z\left[h^{(t-1)},x^{(t)}\right]\right) &\text{update gate}\\ r^{(t)}&= \sigma\left(W_r\left[h^{(t-1)},x^{(t)}\right]\right) &\text{reset gate}\\ \tilde{h}^{(t)}&= \tanh\left( W\left[ r^{(t)}\odot h^{(t-1)},x^{(t)} \right] \right) &\text{candidate state}\\ h^{(t)}&=\left(1-z^{(t)}\right)\odot h^{(t-1)}+z^{(t)}\odot \tilde{h}^{(t)} \\ \end{align*}$

其中$\odot$为elemente-wise product。

update gate $z^{(t)}$会决定最后是保留更多的previous state还是新的candidate state。

reset gate $r^{(t)}$则是决定计算新的candidate state时，参考多少之前的hidden state的信息。

LSTM (Long Short-Term Memory)

更加强力的架构。我们在hidden state $h^{(t)}$的基础上再引入一个新的cell state $c^{(t)}$。

Forget gate:
$f^{(t)} = \sigma\left( W_f\left[h^{(t-1)},x^{(t)}\right] \right)$
决定保留多少之前的cell state信息。
Input gate:
$i^{(t)} = \sigma\left( W_i\left[h^{(t-1)},x^{(t)}\right] \right)$
决定将多少新信息写入新的cell state。
Output gate:
$o^{(t)} = \sigma\left( W_o\left[h^{(t-1)},x^{(t)}\right] \right)$
决定更新hidden state时参考多少cell state的信息。
Update cell state:
$c^{(t)} = f^{(t)}\odot c^{(t-1)} + i^{(t)}\odot \tanh\left( W\left[h^{(t-1)},x^{(t)}\right] \right)$
Update hidden state using cell state:
$h^{(t)} = o^{(t)}\odot \tanh\left(c^{(t)}\right)$

Non-Recurrent Modesl (ConvNets, Transformer)

RNN的结构：

$h^{(t)} = f\left(h^{(t-1)},x^{(t)}\right)$

会带来2个问题：

无法并行计算，因为每个hidden state都依赖前面的state
无法保证长期

但有些序列任务其实不一定需要完整历史，只需要一个固定范围的过去信息。

ConvNets/WaveNet

连续函数上的卷积是：

$(f * g)(t) = \int_{-\infty}^{+\infty} f(\tau)g(t-\tau)\,d\tau$

在图像处理中，我们常见2D（离散）卷积：

$(K * I)(i,j) = \sum_m \sum_n I(i-m,j-n)K(m,n)$

因为sequence是1维的，所以我们可以用1D的ConvNets。

假设有一个长度为 $3$ 的卷积核：

$[k_0,k_1,k_2]$

对于序列：

$x_1,x_2,x_3,x_4,x_5$

普通1D卷积会在序列上滑动窗口。例如某个输出可能是：

$y_t = k_0x_{t-2} + k_1x_{t-1} + k_2x_t$

这和之前的autoregressive model非常像，但ConvNet可以叠很多层、加非线性，从而学习更复杂的模式。

WaveNet是一个用 1D ConvNet 建模 speech / raw audio 的模型。它有两个关键设计：

Causal convolution

确保卷积只能看到之前的信息。
Dilated convolution

（在序列模型里，某个输出 $y_t$ 能看到多少个输入位置，就叫它的receptive field。）

普通卷积看连续相邻输入：
$x_t,x_{t-1},x_{t-2}$
而 dilated convolution 会“隔着看”。

比如 dilation $=2$ 时，可能看：
$x_t,x_{t-2},x_{t-4}$
dilation $=4$ 时，看：
$x_t,x_{t-4},x_{t-8}$
WaveNet常用的做法是每一层dilation翻倍：
$1, 2, 4, 8, \ldots$
这样receptive field会随着层数指数级增长。这样便可以只用少数层就能覆盖很长历史。

Transformer&Attention

Transformer 是非循环模型。它不像RNN那样一步一步更新hidden state，也不像ConvNet那样只看固定窗口，而是用attention让每个token可以直接和其他token交互。

对于当前位置 $x_i$，模型会学习

应该关注序列中哪些位置 $x_j$？
每个位置关注多少？（权重 weight）

（attension中不同深浅的颜色代表不同的重要程度/权重）

Attension有三个核心概念：

$Q,\ K,\ V$

也就是：

Query
Key
Value

可以用一个比喻理解：

Query：我在找什么？
Key：我有什么特征，可以被别人匹配？
Value：如果别人关注我，我实际提供什么信息？

attension的具体计算过程：

第一步，计算 $x_i$ 对每个 $x_j$ 的分数：

$q_i = W_Q x_i$ $\begin{align*} k_j = W_K x_j \end{align*}$ $\mathrm{score}_{ij} = q_i k_j^T$

分数越高，表示$x_i$越应该关注$x_j$。

第二步，对这些分数做 scaling 和 softmax：

$\alpha_{ij} = \operatorname{softmax}_j \left( \frac{q_i k_j^T}{\sqrt{d_k}} \right)$

这里$\alpha_{ij}$是attention weight，表示位置$i$对位置$j$的关注程度。

第三步，用这些权重加权value：

$v_j = W_V x_j$ $z_i = \sum_j \alpha_{ij}v_j$

得到当前位置的新表示$z_i$。

用矩阵的写法则是：

$\begin{align*} Q &= XW_Q \\ K &= XW_K \\ V &= XW_V \\ Z &= \operatorname{softmax} \left( \frac{QK^T}{\sqrt{d_k}} \right) V \\ \end{align*}$

其中$QK^T$是 $n \times n$ 的矩阵。

由于Attention本身无法察觉顺序，所以我们需要加入Positional Encoding，也就是给每个token的embedding加上或拼接一个表示位置的向量。

与RNN相比，Transformer有一个非常重要的特性：训练时可以并行。RNN必须按顺序处理每个token，因为每一步的输出都依赖上一步的隐藏状态；而Transformer没有这种时序依赖，可以一次性把整个序列喂进去同时计算。

但这里有个问题——语言模型在预测某个token时，不应该”提前看到”它后面的内容，否则训练就失去意义了。解决办法是使用causal mask（因果掩码）：在注意力计算中，将每个位置右边的部分遮住，强制模型只能关注当前及之前的token。

以训练句子”I like this movie”为例：

I 只能看 I
like 只能看 I, like
this 只能看 I, like, this
movie 只能看 I, like, this, movie

这四个位置的计算在训练时是同时进行的，效率远高于RNN的逐步处理。

不过，这种并行性只存在于训练阶段。推理时，由于每个新token都依赖前面已生成的结果，模型仍然只能一个token一个token地顺序输出。

Stuctured State Space Models

先来对比一下前面的RNN和Transformer架构：

	Training 训练	Inference 推理
RNN	慢	快（只需要根据当前的hidden state来计算下一个hidden state）
Transformer	快（可并行）	慢（复杂度差不多是序列长度的平方）

那么有没有办法可以设计一个架构同时拥有它们2个的优点呢？

State-space model

考虑输入$u_k \in \mathbb{R}$，输出$y_k \in \mathbb{R}$，以及$\bar{A}\in\mathbb{R}^{N\times N}, \bar{B}\in\mathbb{R}^{N\times 1}, \bar{C}\in\mathbb{R}^{1\times N}$。

$x_k = \bar{A}x_{k-1} + \bar{B}u_k$ $y_k = \bar{C}x_k$

注意到$y_k$可以写成

$\begin{align*} y_k &=\bar{C}\bar{A}^k\bar{B}u_0+\cdots+\bar{C}\bar{A}\bar{B}u_{k-1}+\bar{C}\bar{B}u_k \\ &= K_0 u_0+\cdots + K_{k-1} u_{k-1}+K_k u_k \\ \end{align*}$

等价于$\bar{K}$和$u$的卷积：

$\bar{K} = \left[ \bar{C}\bar{B}, \bar{C}\bar{A}\bar{B}, \bar{C}\bar{A}^2\bar{B}, \ldots, \bar{C}\bar{A}^{L-1}\bar{B} \right]$

重点来了，目前的这个架构有着2种等价的形式，所以同时拥有着2种形式的优势：

递归
$\begin{align*} x_k &= \bar{A}x_{k-1} + \bar{B}u_k \\ y_k &= \bar{C}x_k \end{align*}$
推理时很高效，复杂度为线性。
卷积
$y = \bar{K} * u$
训练时可以并行，因为整段序列的卷积可以一次算。

但是光这样还不够，因为在计算$\bar{K}$时实际上需要计算$\bar{A}^{L-1}$，复杂度依旧很高。

所以我们可以挑选合适的$\bar{A}$的结构来确保复杂度不会太高。

S4 model

S4 (Structed State Space for Sequences) model。

Event Data Modeling

之前考虑的都是数列的值，也就是$x_{t_1},x_{t_2},…$。我们现在来考虑时间（也就是index本身）。

也就是我们现在有一些event，然后我们希望预测下一个event会在什么是发生。

Definition: Temporal Point Processes (TPP) are a class of probabilistic models that describe the distribution of discrete event sequences in continuous time.

TPP defines a generative model for variable-length sequences

$t=\{t_1,\ldots,t_N\}$

on the interval $[0,T]$, with likelihood function $p(\{ t_1,…,t_N \})$.

记历史事件为$\mathcal{H}(t) = \{ t_j < t \} = \{t_j \mid t_j < t \}$。那么模型描述的conditional density（下一个事件发生的时间概率分布）就是

$p^*(t):= p(t \mid \mathcal{H}(t))$

并且不难注意到：

$\Pr\!\left(T_{\text{next event}}\le t\mid \mathcal H(t)\right)=0,$

（严格来讲图里的$\mathcal{H}(t)$应该只画到$t_3$那里。）

CDF：

$F^*(t) = \int_{t_{i-1}}^{t} p^*(u)\,du$

为下一个时间发生在$[t_{i-1},t)$的概率。其中$t_{i-1}$为$t$之前的最后一个事件（的时间）。

Survival function：

$S^*(t) = 1-F^*(t) = \int_t^{\infty} p^*(u)\,du$

为下一个时间不发生在时间$t$之前的概率，或者说下一个时间发生在$t$之后的概率。

conditional intensity function：

$\begin{align*} \lambda^*(t)dt &= \Pr\left(\text{event in }[t,t+dt)\ \mid \text{no event in }[t_{i-1},t)\right) \\ &= \frac{ \Pr\left(\text{event in }[t,t+dt)\ \&\ \text{no event in }[t_{i-1},t)\right) }{ \Pr\left(\text{no event in }[t_{i-1},t)\right) }\\ &= \frac{\int_t^{t+dt} p^*(t)dt}{S^*(t)} \end{align*}$

在已知$[t_{i-1},t)$中没有发生事件的情况下在$[t,t+dt)$中发生事件的概率。

（$\int_t^{t+dt} p^(t)dt$有时候会被省略成$p^(t)dt$。）

也就是：

$\lambda^*(t) = \frac{p^*(t)}{S^*(t)}$

这几个函数之间的关系：

那么该怎么计算一个realization的概率呢？

$\begin{align*} p(\{t_1,t_2,t_3,t_4\}) &= p(t_1)p(t_2 \mid t_1)p(t_3 \mid t_1,t_2)p(t_4 \mid t_1,t_2,t_3)S^*(T) \\ &= p^*(t_1)p^*(t_2)p^*(t_3)p^*(t_4)S^*(T) \\ &= \lambda^*(t_1)\lambda^*(t_2)\lambda^*(t_3)\lambda^*(t_4) \exp\left( -\int_0^T \lambda^*(u)\,du \right) \end{align*}$

Selected TPP Models

**用$\lambda^*(t)$建模**

如果用$\lambda^*(t)$来定义模型，有几个好处：

容易设计具有特定行为的模型，例如全局趋势、爆发性、排斥性。
强度函数比较可解释。
不同过程可以组合，例如两个强度相加。
采样可以更高效。

Homogeneous Poisson Process（HPP）

假设强度恒定：

$\lambda^*(t)=\mu$

服从指数分布：

$p^*(t) = \mu \exp\left( -\int_{t_{i-1}}^{t}\mu\,du \right) = \mu\exp\left(-\mu(t-t_{i-1})\right)$

记$\tau_i = t_i-t_{i-1}$，那么

$P(\tau_i > s) = S(s) = \exp\left(-\mu s \right)$ $\tau_i \sim \text{exponential}(\mu)$

模拟HPP很简单：不断从指数分布采样间隔 $\tau_i$，然后累加得到事件时间：

arrival_times = []
t = 0
while t < T:
    tau ~ Exponential(mu)
    t += tau
    if t < T:
        arrival_times.append(t)

复习一下一个定理：

Theorem （Inverse CDF Transform）:

设

$U \sim \operatorname{Uniform}(0,1)$

设 $F:\mathbb{R}\to[0,1]$ 是某个随机变量的累积分布函数（CDF），定义广义反函数，即quantile function：

$F^{-1}(u):=\inf\{x\in\mathbb{R}:F(x)\ge u\},\quad u\in(0,1).$

令

$X:=F^{-1}(U).$

则 $X$ 的累积分布函数为 $F$，也就是说：

$X\sim F.$

换句话说：

$\Pr(X\le x)=F(x),\quad\forall x\in\mathbb{R}.$

所以取样（Sampling）的话就很简单：

取$u \sim \operatorname{Uniform}(0,1)$，那么

$\tau = F^{-1}(u) = -\frac{1}{\mu}\text{log}(1-\mu)$

HPP 适合描述“事件均匀随机发生”的场景，但表达力很弱。

Inhomogeneous Poisson Process （IPP）

IPP 允许强度随时间变化：

$\lambda^*(t)=g(t)\ge 0$

不依赖历史，但可以捕捉全局趋势。

而事件的数量满足Poisson分布：

$N([a,b]) \sim \operatorname{Poisson} \left( \int_a^b g(t)\,dt \right)$

假如有2个IPP的intensities是$g(t)$和$h(t)$，那么它们的和也是IPP并且intensity就是$g(t) + h(t)$。

Hawkes Process

也叫自激励过程（self-exciting process）。它的intensity等于

$\lambda^*(t) = \mu + \alpha \sum_{t_j\in \mathcal{H}(t)} k_\omega(t-t_j)$

其中

$k_\omega(t-t_i) = \exp\left(-\omega(t-t_i)\right)$

被叫做Triggering kernel。并且参数$\mu, \alpha, \omega \geq 0$。是依赖历史记录的。

Hawkes Process 的核心特点是 self-exciting。一个事件发生后，会暂时提高后续事件的发生概率，然后这个影响逐渐衰减。因此它适合建模 bursty、clustered 的现象，比如社交媒体转发、地震余震、金融市场连锁交易、用户连续发消息等。

参数估计/学习

最大化训练集中所有序列的 log-likelihood：

$\max_{\theta} \sum_{t=\{t_1,\ldots,t_N\}\in \mathcal{D}_{\mathrm{train}}} \log p_{\theta}(\{t_1,\ldots,t_N\})$

其中

$\log p_{\theta}(\{t_1,\ldots,t_N\}) = \sum_{i=1}^{N} \log \lambda^*(t_i) - \int_0^T \lambda^*(u)\,du$

**建模$p^*(t)$**

因为用$\lambda^(t)$建模的话计算log-likelihood部分需要算积分，复杂的情况下会无法计算。所以我们也可以直接用RNN建模$p^(t)$。

利用RNN来记录evnt时间差的历史，并基于这个历史得到$p^*(t)$的分布。

但这样autoregressive TPP还是有几个（和之前RNN类似的）问题：

第一，误差累积。预测下一个事件错了，后面基于它继续预测会越错越远。
第二，计算复杂度高，因为要一步步生成事件。
第三，难以建模长距离依赖，尤其是很早以前的事件对很久之后事件的影响。

再来看个变种：ADD and THIN: Diffusion for Temporal Point Processes。

Noising 阶段：

从真实数据样本 $\mathbf{t}^{(0)}\sim \lambda_0$ 开始，每一步做两件事：

Thin：随机删除一部分已有事件。
Add：从 HPP 添加一些噪声事件。

经过很多步后，原始事件序列逐渐变成一个简单的 HPP 噪声分布：

$\mathbf{t}^{(N)}\sim \lambda_{HPP}$

Denoising 阶段：

从噪声事件序列开始，逐步采样回更干净的事件序列，最终恢复数据分布：

$\mathbf{t}^{(0)}\sim \lambda_0$

Graphs

Ranking

核心问题：该怎么给节点排序？

最典型的例子就是 Web 网页排序。互联网可以看成一个巨大的有向图：网页是节点，超链接是边。问题是，面对海量网页，怎么判断哪个网页更“重要”？

PageRank

PageRank 的核心递归思想是：A page is important if many importants pages point to it。如果有很多重要页面指向某个页面，那么这个页面就很重要。

这是一种“投票”机制：每个页面会把自己的重要性分给它指向的页面；如果一个页面有很多出边，它的投票会被分摊。也就是说，如果页面 $i$ 的 PageRank 是 $r_i$，出度是 $d_i$，并且它指向页面 $j$，那么它给 $j$ 的贡献是：

$\frac{r_i}{d_i}$

因此页面 $j$ 的排名分数是所有指向它的页面贡献之和：

$r_j = \sum_{i \to j} \frac{r_i}{d_i}$

这里 $d_i$ 表示节点 $i$ 的出度。

直观理解：
如果一个非常重要的网站只链接到少数几个页面，那么它给每个被链接页面的“投票”很强；如果一个页面本身不重要，或者它链接出去很多页面，那么它对单个页面的贡献就小。

不能发现最后所有的等式会组合成一个方程组。为了更好的表达，我们引入随机邻接矩阵M：

$M_{ji} = \begin{cases} \dfrac{1}{d_i}, & \text{if } i \to j,\\ 0, & \text{otherwise}. \end{cases}$

而原来的 PageRank 方程可以写成：

$r = Mr$

也就是说，PageRank向量$r$是矩阵$M$的一个特征向量，对应特征值为 $1$。

因为M是列随机矩阵（每一列的和为 $1$），所以最大的特征值就是1。那么现在问题变成了求M的最大特征值的特征向量。

可以用幂迭代：

初始化一个归一化的随机向量 $r$；
不断更新：
$r \leftarrow M r$
直到收敛。

并且M是稀疏矩阵（sparse），实际计算时只需要考虑每个节点的入邻居，而不是构造一个完整的巨大矩阵。

还可以把网页跳转理解成是一个markov chain，下一步跳转的网站只取决于当前在哪个网页：

$\text{Pr}(X_{t+1}=j \mid X_t = i) = M_{ji}$

为了保证随机游走有唯一稳定的长期分布，图需要满足2个条件：

不可约，irreducible：从任意状态都可以到达任意其他状态。
换成网页图，就是从任何网页出发，都有路径能到达其他网页。
非周期，aperiodic：状态不会只能以固定周期返回自己。
如果一个状态只能每隔 3 步、6 步、9 步回到自己，那么它就是周期性的。

在irreducible的Markov chain里，只要有一个状态是非周期的，那么所有状态都是非周期的。

普通的PageRank 会遇到的三个问题。

第一，dead ends，死胡同节点。有些网页没有出链，随机游走走到这里就“无路可走”，导致重要性泄漏出去。

第二，spider traps，蜘蛛陷阱。某一组网页的所有出链都指向组内，随机游走一旦进入这个集合，就出不来了，最终这些节点会吸收所有重要性。

第三，periodic states，周期状态。如果从某个状态出发只能以固定周期回到它自己，随机游走可能不会稳定收敛，而是来回震荡。

所以，单纯沿链接随机游走还不够，需要一个修正机制。

Random Teleports

随机冲浪者在每一步有两个选择：

$\beta$的概率沿着一个随机出链继续走；
$1-\beta$的概率跳到某个随机页面。

于是 PageRank 公式变为：

$r_j = \sum_{i \to j} \beta \frac{r_i}{d_i} + (1-\beta)\frac{1}{N}$

矩阵形式为：

$A = \beta M + (1-\beta)\frac{1}{N}\mathbf{1}\mathbf{1}^T$

最终求：

$r = A r$

带随机传送后，矩阵 $A$ 会变成稠密矩阵，因为随机传送意味着每个节点都可能跳到任意节点。所以不能用稠密矩阵来直接来计算PageRank，而是先计算稀疏M的部分，然后再加上常数项。

但PageRank也有一些问题：

PageRank 衡量的是通用流行度，不一定适合特定主题。例如一个体育主题下的权威网页，未必在全网 PageRank 很高。

解决思路是Topic-Sensitive PageRank。
PageRank 容易受到 link spam 影响。有人可以人为制造链接结构，提高目标网页 PageRank。

解决思路是TrustRank。
PageRank 只有单一重要性分数，但现实中“重要性”可能有多种定义。
解决思路包括Hubs-and-Authorities (HITS, Hyperlink-Induced Topic Search)等。

Topic-Sensitive PageRank的思想是：随机传送时，不是跳到所有页面，而是跳到某个主题相关页面集合 $S$。例如，如果用户搜索体育内容，就让随机冲浪者更可能跳到体育相关页面；如果用户搜索历史内容，就跳到历史相关页面。这样得到的 PageRank 不再是全局流行度，而是主题相关重要性。

矩阵形式可以写成：

$r = \beta M r + (1-\beta)\pi$

其中 $\pi$ 是 teleport vector。如果是主题集合 $S$，那么：

$\pi_i = \begin{cases} \frac{1}{|S|}, & i \in S \\ 0, & otherwise \end{cases}$

还可以进一步推广到任意 teleport vector $\pi$，只要满足 $\sum_i \pi_i = 1$即可。

TrustRank的思想是：用可信网页作为 teleport set，做一种 topic-sensitive PageRank。例如使用 .edu 域名，或者非美国学校的类似可信域名。这样随机传送主要从可信页面出发，垃圾链接农场就不容易凭空获得很高权重。

Clustering

输入是一个图 $G=(V,E)$，目标是把顶点 $V$ 分成若干个 cluster。图聚类主要根据对象之间的关系聚类，也就是边所表达的连接信息。

一个“好”的图聚类通常满足两点：同一个簇内部有很多边，也就是 intra-cluster edges 多；不同簇之间只有少量边，也就是 inter-cluster edges 少。

Clustering的算法大致可以分成2种：

partitioning approaches 要求每个节点恰好属于一个簇；
non-partitioning approaches 允许簇重叠，也允许有离群点不属于任何簇。

我们这里主要关注weighted undirected graphs 上的 non-overlapping clustering。

可以把问题转换成优化问题：给定图 $G=(V,E)$，把顶点集 $V$ 分成 $k$ 个簇 $\mathcal C=\{C_1,\dots,C_k\}$，

$\text{min}Q(\mathcal C) \text{ s.t } \biguplus_{i=1,...,k}C_i = V$

（disjoint union）。

Minimum Cut

考虑2个cluster的情况：把图分成两部分 $C_1,C_2$，让跨越两个簇之间的边权总和尽可能小。这个量叫 cut：

$\text{cut}(C_1,C_2)=\sum_{v_1\in C_1, v_2\in C_2} w(v_1,v_2)$

但是minimum cut 有严重缺陷：它倾向于把很小的一组节点，甚至单个节点，从大图中切出来，因为这样跨簇边可能最少；而且它只关心 inter-cluster edges，没有考虑簇内部是否紧密。

所以可以引入归一化cut标准：

$\begin{align*} \text{RatioCut}(C_1,C_2) &= \frac{\text{cut}(C_1,C_2)}{|C_1|} + \frac{\text{cut}(C_2,C_1)}{|C_2|} \\ \text{NormalizedCut}(C_1,C_2) &= \frac{\text{cut}(C_1,C_2)}{\text{vol}(C_1)} + \frac{\text{cut}(C_2,C_1)}{\text{vol}(C_2)} \end{align*}$

其中 $\text{vol}(C)$ 是节点集合的 volume，即集合内所有节点度数之和。

但是当cluster数量上来了之后这个问题便变成了NP-hard的问题。

换个思路，我们可以把它变成一个labeling的问题：对每个节点 $v$，赋予一个标签 $f_v$ 表示它属于哪个簇。二分类时，可以令 $f_v\in\{+1,-1\}$。聚类目标是让标签函数 $f$ 在图上“平滑”：相邻节点的标签应该相似。于是可以最小化：

$\sum_{\{u,v\}\in E} W_{uv}(f_u-f_v)^2$

如果一条边权重大，而两个端点标签差异大，那么惩罚就大。这个目标强迫连接紧密的点有相似标签。

定义图拉普拉斯矩阵（Laplacian matrix）：

$L=D-W$

即

$L_{uv} = \begin{cases} -W_{uv}, & \text{if } (u,v) \in E \\ \deg(v), & \text{if } u = v \\ 0, & \text{otherwise} \end{cases}$

其中 $W$ 是加权邻接矩阵，$D$ 是（加权）度矩阵（$D_{ii}=\text{deg}(i)=\sum_jW_{ij}$）。对无向图，$W$ 对称，$D$ 是对角矩阵。

拉普拉斯矩阵（Laplacian matrix）有几个重要性质：

$f^\top L f = \sum_{\{u,v\}\in E} W_{uv}(f_u-f_v)^2$

刚好可以把我们刚才的目标函数（objective function）转改写成矩阵的乘法。

证明：

注意到：

$\begin{align*} f^\top Df = \sum_u D_{uu}f_u^2 &= \sum_u \deg(u)f_u^2 \\ &= \sum_u \left(\sum_v W_{uv}\right)f_u^2 \\ &= \sum_{u,v} W_{uv}f_u^2 \end{align*}$

所以

$\begin{align*} f^\top Lf &= f^\top Df - f^\top Wf\\ &= \sum_{u,v} W_{uv} f_u^2 - \sum_{u,v} W_{uv} f_u f_v \\ &= \frac12 \left( \sum_{u,v} W_{uv}f_u^2 + \sum_{u,v} W_{uv}f_v^2 - 2\sum_{u,v} W_{uv}f_uf_v \right) \\ &= \frac12 \sum_{u,v} W_{uv} \left( f_u^2-2f_uf_v+f_v^2 \right) \\ &= \frac12 \sum_{u,v} W_{uv}(f_u-f_v)^2 \end{align*}$

$\square$

2. Laplacian是additive的：

$L_{G \cup H} = L_G + L_H$

3. $L$ 是对称半正定矩阵（symmetric, positive semi-definite）

4. $L$ 的最小特征值是 0，对应的特征向量是全 1 向量 $\mathbf 1$，或者任意常数倍 $c\mathbf 1$。

证明0是最小的特征值：

设 $x$ 是 $L$ 的任意一个特征向量，对应特征值是 $\lambda$：

$Lx=\lambda x$

两边左乘 $x^\top$：

$x^\top Lx = x^\top \lambda x$

也就是：

$x^\top Lx = \lambda x^\top x$

所以：

$\lambda = \frac{x^\top Lx}{x^\top x} \geq 0$

5. $L$ 有 $n$ 个非负实特征值：

$0=\lambda_1\le \lambda_2\le \cdots \le \lambda_n$

来看2个cluster的情况。考虑min RatioCut：

$\min_{C_1 \subset V} \frac{\operatorname{cut}(C_1,C_2)}{|C_1|} + \frac{\operatorname{cut}(C_2,C_1)}{|C_2|}, \qquad \text{where } C_2 = V \setminus C_1 := \bar{C}_1$

定义一个比较奇怪的indicator vector $f_{C_1}$：

$f_{C_1}[i] = \begin{cases} +\sqrt{\frac{|\bar{C}_1|}{|C_1|}}, & \text{if } v_i \in C_1 \\ -\sqrt{\frac{|C_1|}{|\bar{C}_1|}}, & \text{else} \end{cases}$

但它有些非常重要的性质：

1. $\sum_i f_{C_1}[i] = 0$ ，也就是$f_{C_1} \perp \vec{1}$

2. $f_{C_1}^T \cdot f_{C_1} = \|f_{C_1}\|_2^2 = |V|$ （长度是常量）

$f_{C_1}^T \cdot L \cdot f_{C_1} = \cdots = |V| \cdot \left[\frac{\operatorname{cut}(C_1,C_2)}{|C_1|} + \frac{\operatorname{cut}(C_1,C_2)}{|\bar{C}_1|}\right] = |V| \cdot \text{ratio cut}$

当我们给$f$的要求放宽松一点，只要求它满足1，2这两条性质，那么就可以把问题变成：

$\min_{f_{C_1} \in \mathbb{R}^{|V|}} \left\{ f_{C_1}^T \cdot L \cdot f_{C_1} \right\} \quad \text{subject to } f_{C_1} \perp \vec{1} \text{ and } \|f_{C_1}\|_2 = \sqrt{|V|}$

注意到，对于所以的对称矩阵$L$都有：

$\begin{align*} \lambda_1 &= \min_{\|x_1\|=1} x_1^T L x_1 \\ \lambda_2 &= \min_{\substack{\|x_2\|=1, \; x_2 \perp x_1}} x_2^T L x_2 \\ \lambda_3 &= \cdots \end{align*}$

所以上面优化问题实际上要找到$f_{C_1}$实际上就是$L$第二小的eigenvector。

那么最后怎么通过$f_{C_1}$来还原cluster呢？

2个cluster的情况直接考虑$f_{C_1}$的符号即可
多个cluster就需要基于$f_{C_1}$的值做一个k-means clustering。

现在再推广到cluster数量$k>2$的情况。

对每个簇 $C$，定义一个 indicator vector：

$h_C(i)= \begin{cases} \frac{1}{\sqrt{|C|}}, & v_i\in C\\[4pt] 0, & \text{otherwise} \end{cases}$

然后把这些向量作为列拼成矩阵：

$H= [h_{C_1},h_{C_2},\dots,h_{C_k}]$

不难注意到：

$H^TH=I$ $(H^T L H)_{ii} = h_{C_i}^T \cdot L \cdot h_{C_i} = \frac{\operatorname{cut}(C_i, V \setminus C_i)}{|C_i|}$

所以

$\operatorname{RatioCut}(C_1, \ldots, C_k) = \sum_{i=1}^k \frac{\operatorname{cut}(C_i, V \setminus C_i)}{|C_i|} = \sum_{i=1}^k (H^T L H)_{ii} = \operatorname{trace}(H^T L H)$

同样，问题就变成了：

$\min_{C_1,\ldots,C_k} \operatorname{trace}(H^T L H) \quad \text{subject to } H^T H = \mathrm{Id} \text{ and } H \text{ as defined before}$

如果我们给H的条件放宽松一点，那么就是

$\min_{H \in \mathbb{R}^{|V| \times k}} \operatorname{trace}(H^T L H) \quad \text{subject to } H^T H = \mathrm{Id}$

类似的，最优解H就是L的前k个最小eigenvectors，因为：

设 $L$ 的特征分解为：

$L=U\Lambda U^\top$

其中：

$\Lambda= \text{diag}(\lambda_1,\dots,\lambda_n)$

并且：

$\lambda_1\le \lambda_2\le\cdots\le\lambda_n$

目标函数是：

$\text{trace}(H^\top LH)$

代入特征分解：

$\text{trace}(H^\top U\Lambda U^\top H)$

令：

$Y=U^\top H$

因为 $U$ 是正交矩阵，且 $H^\top H=I$，所以：

$Y^\top Y=I$

目标变成：

$\text{trace}(Y^\top \Lambda Y) = \sum_{j=1}^k y_j^\top \Lambda y_j = \sum_{j=1}^k \lambda_1 y_{j1}^2+ \cdots \lambda_n y_{jn}^2$

并且

$\forall j: \|y_j\|^2=1$

所以$\text{trace}(Y^\top \Lambda Y)$的最优解就是

$y_1 = e1,...,y_k = e_k$

然后再带回去原问题：

$H^\star = U[e_1,e_2,\dots,e_k] = [u_1,u_2,\dots,u_k]$

spectral embeddings

$H \in \mathbb{R}^{n \times k}$的列向量是L的前k个eigenvector，然后行向量就是每个node的embedding。

按照一开始严格的H要求，得出来的每个行向量应该只有1项是不为0的，表示cluster。即便是后面没有严格要求了，这些向量依旧能携带cluster信息。

可以针对这些embedding做k-means确定每个的label/cluster。

Modularity

前面cut-based 方法的一个局限是：通常需要提前知道簇数 $k$。

而modularity 的核心思想是：不直接看簇内边有多少，而是比较“实际簇内边数”和“在随机图里期望的簇内边数”差多少。

Probabilistic Approaches

Prediction

图上的预测问题可以分成几个层级：

节点级任务（node-level）：预测单个节点的属性或标签。
图级任务（graph-level）：预测整张图的属性。
边级任务（edge-level）：可以预测已有边的属性，也可以预测新边是否可能存在。

学习任务也可以分成几种：

Transductive learning：只预测当前给定的未标注样本

给定一批有标签训练样本 $T$，一批无标签测试样本 $U$，以及可能存在的额外知识，比如图结构 $W$。目标不是学习一个可以泛化到任意新样本的函数，而是只预测当前这批无标签样本 $U$ 的标签。
Inductive learning：学习一个可泛化的预测函数

学习一个映射 $f:\mathcal X \to \mathcal Y$，它可以应用到未来新来的样本 $x_{new}$。图级任务往往是 inductive 的。

此外还有一个概念Semi-supervised learning：同时利用有标签和无标签数据训练。

transductive learning 几乎总是 semi-supervised，因为我们明确拿到了有标签节点和无标签节点，并且目标就是利用所有这些数据来预测无标签节点。

但图是离散结构，传统机器学习通常处理固定长度向量。那么，怎样把图上的节点、边、整张图变成机器学习模型能处理的形式？

答案主要有两类：
第一类是 embedding，显式地把节点或图映射成向量；第二类是 kernel methods，不一定显式生成向量，而是直接定义两个对象之间的相似度函数。

Embedding

Node embedding：把图变换成一种表示，使每个节点都对应一个向量。形式上，节点嵌入函数是

$\Phi: V \to \mathbb R^d$

也就是把每个节点映射到 $d$ 维实数空间。embedding 空间里距离近的节点，应当在图结构意义上“相似”。

embedding 没有唯一标准。不同方法的差异在于：如何定义 embedding 空间中的“相似”，以及想捕捉图的哪种性质。

可以粗略地分为2类：

Spectral embeddings：节点相似通常意味着它们属于相同 cluster。
Deep embeddings：相似性取决于 loss function 和模型结构。

Spectral embeddings

详见前面clustering那章的Spectral embeddings。

Deep embeddings

DeepWalk：把图变成一组 random walks（也就是sequence of tokens），然后在这些随机游走序列上训练 word2vec。对每个节点 $v_i$，采样多条随机游走；这些 random walks 相当于 word2vec 里的句子，而节点 ID 相当于词。

训练完成后，在图中距离近、共享邻居多的节点，会得到相似的向量表示。

Graph2Gauss：把节点映射成embedding空间中的一个高斯分布。目标是让一个节点的一阶邻居比二阶邻居更近，二阶邻居又比更远的节点更近，也就是保留图中的距离排序。（暂时先不考虑边权）

它学习一个映射：

$f_\theta: v_i \mapsto \mathcal N(\mu_\theta(x_i), \Sigma_\theta(x_i))$

对于 attributed graph，$x_i$ 可以是节点特征；对于 plain graph，$x_i$ 可以是节点 ID，比如 one-hot 编码。

mean $\mu_i$表示节点在embedding space中的大概位置，而variance $\Sigma_i$表示不确定性。

用 KL divergence 来衡量两个节点对应的高斯分布的差别有多大（embedding的距离）：

$E_{uv} = KL(f_\theta(u) \| f_\theta(v))$

Graph2Gauss 的 inductive 能力依赖于 $x_i$。

如果 $x_i$ 是真实节点属性，那么模型可以泛化到新节点。
但如果 $x_i$ 只有 node ID 或 one-hot encoding，那么模型更多是在记住训练图里的节点，泛化能力就弱很多。

更进一步，如果两个节点 $u$ 和 $v$ 有完全相同的属性 $x$，但在图中距离很远，那么会出现一个问题：

如果输入只看属性 $x$，它们输入一样；
同一个函数 $f_\theta$ 对一样的输入可能输出一样或非常相近的高斯分布；
但图结构却要求它们离得远。

所以核心trade-off是：

用节点属性做输入，inductive能力强，但可能区分不了属性相同、结构位置不同的节点；
用node ID / one-hot做输入，可以区分节点，但inductive泛化能力变弱。

Graph embedding

自然想法是：先算出每个节点的 embedding，再把所有节点 embedding 聚合成一个 graph embedding。问题是，不同图的节点数不同，直接拼接节点向量会导致维度不一致。PPT 给出的解决方法是：聚合所有节点 embedding，例如 mean pooling。

也就是说，如果一张图有 $n$ 个节点，每个节点embedding是 $d$ 维，那么可以取平均：

$h_G = \frac{1}{n}\sum_{i=1}^{n} h_i$

其中

$h_i \in \mathbb R^d$

表示图$G$中第$i$个节点的embedding；

这样无论图有多少节点，最终都得到一个固定长度的 $d$ 维图向量。然后就可以把这个 graph embedding 输入 SVM、神经网络等模型。

除此之外还有其他聚合的方法：

$h_G^{sum} = \sum_{i=1}^{n}h_i$

保留规模信息，因为节点越多，sum 通常越大。

$h_G^{max} = \max_{i=1}^{n} h_i$

保留某些最强激活特征，但也不直接保留数量。

Kernel methods

很多算法只依赖内积 $\phi(x)^T\phi(x’)$，所以我们不需要真的算出高维特征 $\phi(x)$，只需要找到一个 $k(x,x’)$ 满足：

$k(x,x')=\phi(x)^T\phi(x')$

然后后续直接计算$k(x,x’)$即可。不再需要先embedding然后再计算inner product。

Graph kernel methods的思路非常类似：与其像 embedding 方法那样显式学习 $\phi(G)$，不如直接定义两个图之间的相似度函数：

$k: \mathcal G \times \mathcal G \to \mathbb R$

也就是输入两张图，输出它们有多相似。

graph kernels 通常通过共享结构模式来衡量相似性，例如 random walks、substructures、Weisfeiler-Lehman labels。它们主要用于 graph-level tasks，但也存在 node-level kernel。

Label Propagation

之前是先把图转换成向量，然后再用机器学习的模型。我们也可以直接在图结构上做节点分类/transductive node classification：给定一张图，其中一部分节点有标签，另一部分节点没有标签，目标是在同一张图上预测未标记节点的标签。

核心假设叫homophily/assortativity，也就是物以类聚，人以群分”：如果两个节点有边相连，那么它们更可能拥有相同标签。

先考虑二分类情况。formal的定义：

- 节点$V = S \cup U$（S: seed，labeled；U: unlabeled）

- 对称加权邻接矩阵（Symmetric weighted adjacency matrix）$W \in \mathbb{R}^{|V| \times |V|}$

- 已标记节点 $i \in S$ 的标签记作 $\hat y_i \in \{0,1\}$

- 希望预测所有节点的标签 $y_i$，$i \in V$。

和之前cluster的思路类似，定义Energy：

$E(\mathbf y)=\frac{1}{2}\sum_{ij} w_{ij}(y_i-y_j)^2$

那么整体的问题便变成了一个constrained integer optimization problem：

$\min_{\mathbf y\in \{0,1\}^{|V|}} \frac{1}{2}\sum_{ij}w_{ij}(y_i-y_j)^2 \quad \text{s.t. } y_i=\hat y_i,\ \forall i\in S$

也可以改写成图拉普拉斯矩阵的形式：

$\min_{\mathbf y\in \{0,1\}^{|V|}} \mathbf y^T L\mathbf y \quad \text{s.t. } y_i=\hat y_i,\ \forall i\in S$

其中：

$L = D-W$

是graph Laplacian，$D$ 是度矩阵，$W$ 是邻接矩阵。

当把图拉普拉斯矩阵按照 labeled nodes 和 unlabeled nodes 分块：

$L= \begin{bmatrix} L_{SS} & L_{SU}\\ L_{US} & L_{UU} \end{bmatrix}$

标签向量也相应分块：

$\mathbf y = \begin{bmatrix} \mathbf y_S\\ \mathbf y_U \end{bmatrix} = \begin{bmatrix} \hat{\mathbf y}_S\\ \mathbf y_U \end{bmatrix}$

其中 $\hat{\mathbf y}_S$ 是已知的 seed labels，$\mathbf y_U$ 是要学习的未标记节点标签。

在$L_{UU}$可逆的情况下，不难得到最优解：

$\mathbf y_U = - L_{UU}^{-1}L_{US}\hat{\mathbf y}_S$

证明：

首先

$\mathbf y^T L\mathbf y = \mathbf y_S^T L_{SS}\mathbf y_S + \mathbf y_S^T L_{SU}\mathbf y_U + \mathbf y_U^T L_{US}\mathbf y_S + \mathbf y_U^T L_{UU}\mathbf y_U$

因为 $L$ 是对称矩阵，所以

$\mathbf y^T L\mathbf y = \mathbf y_S^T L_{SS}\mathbf y_S + 2\mathbf y_U^T L_{US}\mathbf y_S + \mathbf y_U^T L_{UU}\mathbf y_U$

因为$\mathbf y_S$是已知的，所以其中真正需要最小化只有：

$\mathbf y_U^T L_{UU}\mathbf y_U + 2\mathbf y_U^T L_{US}\hat{\mathbf y}_S$

对$\mathbf y_U$求偏导得到：

$\begin{align*} \nabla_{\mathbf y_U} \left( \mathbf y_U^T L_{UU}\mathbf y_U + 2\mathbf y_U^T L_{US}\hat{\mathbf y}_S \right) &= 0 \\ 2L_{UU}\mathbf y_U + 2L_{US}\hat{\mathbf y}_S &=0\\ \mathbf y_U &= - L_{UU}^{-1}L_{US}\hat{\mathbf y}_S \end{align*}$

$\square$

而想要$L_{UU}$可逆也并不困难：

Theorem:

$L_{UU}\text{ 可逆}\iff\text{每个连通分量都至少包含一个 labeled node}$

证明：

对任意 $\mathbf z\in\mathbb R^{|U|}$，有：

$\mathbf z^T L_{UU}\mathbf z = \begin{bmatrix} 0\\ \mathbf x \end{bmatrix}^T L \begin{bmatrix} 0\\ \mathbf x \end{bmatrix}$

W.l.o.g这里可以把 labeled nodes 上的值设成 0，unlabeled nodes 上的值设成 $\mathbf x$。

由于 graph Laplacian 的二次型是边上差异的平方和：

$\mathbf z^T L\mathbf z = \frac{1}{2}\sum_{i,j}w_{ij}(z_i-z_j)^2$

所以：

$\mathbf x^T L_{UU}\mathbf x = \frac{1}{2}\sum_{i,j}w_{ij}(z_i-z_j)^2$

其中：

$z_i= \begin{cases} 0, & i\in S\\ x_i, & i\in U \end{cases}$

考虑$\mathbf x^T L_{UU}\mathbf x=0$的情况：

$\frac12\sum_{i,j}w_{ij}(z_i-z_j)^2=0$

意味着每一条正权重边两端的 $z$ 值都必须相等。

如果每个连通分量都至少包含一个labeled node。那么这个连通分量里的所有node的label也都需要等于这个 labeled node，也就是$x=z=0$。所以$L_{UU}$可逆。

如果某个连通分量没有 labeled node，那么就存在一个$x \neq 0$使得$\mathbf x^T L_{UU}\mathbf x$，代表$L_{UU}$不可逆。

$\square$

那有$K$个类别的情况该怎么办？答案是使用 one-hot notation：

$y_{ik}= \begin{cases} 1, & \text{if node } i \text{ is of class } k\\ 0, & \text{else} \end{cases}$

也就是说，每个节点不再只有一个标量标签，而是一个 $K$ 维标签向量。多分类下的能量函数变成：

$E(Y)=\sum_{i,j} w_{ij}(\mathbf y_i-\mathbf y_j)^T(\mathbf y_i-\mathbf y_j)$

经典 Label Propagation 通常假设 homophily，也就是相连节点倾向于同类。但图中也可能存在其他网络效应，例如 heterophily，即“异类相吸”。

为此，可以引入一个compatibility matrix：

$H$

然后能量函数变为：

$E(Y)=\sum_{i,j} w_{ij} (\mathbf y_i-\mathbf y_j)^T H (\mathbf y_i-\mathbf y_j)$

Soft constraint

我们前面考虑的都是hard constraint：seed 节点标签必须完全等于给定标签。

但我们还可以考虑soft constraint的版本，即seed labels 不再是绝对固定，而是通过一个正则项鼓励模型接近它们：

$\min_{F \in \mathbb{R}^{N \times K}} \left\{ \sum_{i=1}^{N} \sum_{j=1}^{N} W_{i,j} \left\| d_i^{-1} F_{i*} - d_j^{-1} F_{j*} \right\|^2 + \mu \sum_{i=1}^{N} d_i^{-1} \left\| F_{i*} - Y_{i*} \right\|^2 \right\}$

目标函数里有两部分：

第一部分仍然是图平滑项，要求连接节点的预测结果接近；

第二部分是 label fitting 项，要求预测结果 $F$ 接近给定标签矩阵 $Y$。

Y每行是 seed nodes 的 one-hot 向量，非 seed 节点对应零向量；$\mu$ 是正则化参数；$d_i$ 是节点度。

其解析解是：

$F=(1-\alpha)(I_N-\alpha WD^{-1})^{-1}Y$

其中：

$\alpha=\frac{2}{2+\mu}$

Attributed Graphs

前面的Label Propagation主要依赖图结构和部分标签。但现实中很多图还有节点特征，例如社交网络用户的年龄、文本、兴趣向量，论文网络中文章的关键词向量等。

这种图为attributed graphs：每个节点除了图结构之外，还有一个特征向量：

$X_v\in \mathbb R^d$

同时部分节点有标签 $Y$。问题变成：能不能同时利用节点特征、标签和图结构？

PPNP

PPNP的全称是Personalized Propagation of Neural Predictions。

它的核心思想是：

Predict then Propagate，也就是先用节点特征做预测，再把预测结果沿图扩散。

第一步：Predict。

对每个节点 $i$，用神经网络 $f_\theta$ 根据节点特征 $X_i$ 预测初始类别概率：

$H_i^{(0)}=f_\theta(X_i)$

这里 $f_\theta$ 是可学习的 neural network。

第二步：Propagate。

不是传播原始标签，而是传播神经网络的预测结果 $H^{(0)}$。最后的类别预测结果：

$Z= \text{softmax} \left( \alpha \left( I_n-(1-\alpha)\tilde D^{-1/2}\tilde A\tilde D^{-1/2} \right)^{-1} H^{(0)} \right)$

其中 $\tilde A$ 和 $\tilde D$ 是加入 self-loops 后的邻接矩阵和度矩阵。

注意到这部分可以展开成无穷级数：

$S:=\tilde D^{-1/2}\tilde A\tilde D^{-1/2}$ $\left(I-(1-\alpha)S\right)^{-1} = I+(1-\alpha)S+(1-\alpha)^2S^2+(1-\alpha)^3S^3+\cdots$

直观上，PPNP做了两件事：

神经网络负责看节点自己的特征；
Personalized PageRank扩散负责结合图结构，让相邻或相近节点的预测互相影响。

PPNP的训练目标是：只在有标签的节点上计算 $Z$ 和真实标签 $Y$ 之间的 cross entropy loss。然后通过backpropagation端到端训练神经网络参数 $\theta$。

APPNP

Approximate Personalized Propagation of Neural Predictions。

PPNP 的精确解里有矩阵逆：

$\left( I_n-(1-\alpha)\tilde D^{-1/2}\tilde A\tilde D^{-1/2} \right)^{-1}$

大图上直接算这个很贵。所以可以利用随机游走迭代近似精确 Personalized PageRank 解：

$H^{(l+1)} = (1-\alpha)\tilde D^{-1/2}\tilde A\tilde D^{-1/2}H^{(l)} + \alpha H^{(0)}$

GNN

很多图预测任务不能只看节点特征，还必须看复杂的结构关系。

GNN的核心idea：让模型端到端地学习节点特征和图结构应该如何结合。

一个核心约束是Permutation Invariance，置换不变性。同一张图如果只是重新编号节点，本质上没有变；因此图级预测模型的输出不应该依赖任意节点排序，即应满足：

$M(G)=M(\pi(G))$

不能简单把邻接矩阵当图片处理，因为节点顺序是任意的，而且也不希望引入 $O(n^2)$ 的复杂度。

CNN中的idea是：在图像中，一个像素的新表示受自己和邻居像素影响，同时卷积核参数共享。

可以把它泛化到图卷积NN（Convolutional Neural Networks for Graphs）：一个节点的新表示应该由它自己和邻居的旧表示共同决定，并且所有节点共享同一套更新规则。

和sequences/images不同，图上可以定义很多种不同的卷积。目前文献中最常见两类是：

Spatial GNNs：空间视角，直接在图上做可微消息传递。
Spectral GNNs：谱视角，借助图信号处理和谱滤波定义卷积。

Spatial GNNs

Spatial GNN 的核心是 message passing，消息传递。对于图 $G=(V,E,X)$，每个节点 $v$ 在第 $t$ 层有隐藏表示 $h_v^{(t)}$，初始表示通常是节点特征：

$h_v^{(0)}=X_v$

每一层更新为：

$h_v^{(t)} = \text{upd}^{(t)} \left( h_v^{(t-1)}, \text{AGG}\{h_w^{(t-1)} \mid w \in N(v)\} \right)$

节点$v$收集邻居$w$的上一层表示，先聚合，再和自己的旧表示一起更新。

$upd^{(t)}$ 可以是任意可微函数，例如神经网络
同一层中所有节点共享权重
这个过程不依赖节点顺序

第$K$层之后，节点表示会包含$K$-hop邻域的信息。

所以GNN的深度（depth）同时有2种理解：

transformation depth，即非线性变换堆叠多少层
propagation depth，即信息在图上传播了多远

2个实例：

1. GCN 使用带度归一化的邻居求和：

$h_v^{(t)} = \text{ReLU} \left( W^{(t)} \sum_{w\in N(v)\cup\{v\}} \frac{h_w^{(t-1)}}{\sqrt{(\deg(w)+1)(\deg(v)+1)}} \right)$

度归一化的目的是防止高度节点的聚合值过大，也和归一化 Laplacian 有关。

可以写成矩阵形式：

$H^{(t)} = \sigma \left( \tilde{D}^{-1/2} \tilde{A} \tilde{D}^{-1/2} H^{(t-1)} W^{(t-1)} \right)$

其中 $\tilde{A}=A+I_n$，即加了自环；归一化传播矩阵可以避免梯度爆炸。

2. GraphSAGE 则是把“自身信息”和“邻居均值信息”分别做变换：

$h_v^{(t)} = \text{ReLU} \left( W_1^{(t)}h_v^{(t-1)} + W_2^{(t)} \cdot \text{mean}_{w\in N(v)} h_w^{(t-1)} \right)$

它的特点是把 self-message 和 neighbor-message 分开建模。

Spectral GNNs

先定义归一化 Laplacian：

$\hat{L}=I-D^{-1/2}AD^{-1/2}$

并进行特征分解：

$\hat{L}=U\Lambda U^T$

其中 $U$ 的列是 Laplacian 的特征向量，$\Lambda$ 是特征值对角矩阵。

图信号可以看作定义在节点上的函数 $x:V\to \mathbb{R}$，即一个向量 $x\in \mathbb{R}^N$。Graph Fourier transform (GFT)定义为：

$\tilde{x}=U^Tx$

逆变换（inverse GFT）为：

$x=U\tilde{x}$

图上的谱滤波（spectral filtering on graphs）则写作：

$y=Ug(\Lambda)U^Tx$

其中$g$是spectral filter，会把$\Lambda$映射到对角矩阵。

信号处理中的傅里叶的那些$cos(k\pi)$和$sin(k\pi)$是orthonormal basis，矩阵$U$里的列向量也是orthonormal basis。

Spectral GNN的main idea就是把spectral filter$g_\theta$当作参数。每层的更新公式：

$H^{(l+1)} = \phi \left( Ug_\theta(\Lambda)U^T \varphi(H^{(l)}) \right)$

Spectral propagation$Ug_\theta(\Lambda)U^T$ 可以一次覆盖多个 hop，因此和逐层局部传播的message passing视角不同。

但直接做特征分解代价高，所以常用多项式滤波器：

$g_\theta(\Lambda) = \sum_{k=0}^K \theta_k \Lambda^k$

这样可以避免直接使用U来计算：

$H^{(l+1)} = \phi\left(g_\theta(\hat{L})\varphi\left(H^{(l)}\right)\right)$

因为

$U g_\theta(\Lambda) U^T = \sum_{k=0}^{K} \theta_k U \Lambda^k U^T = \sum_{k=0}^{K} \theta_k (U \Lambda U^T)^k = \sum_{k=0}^{K} \theta_k \hat{L}^k$

Spectral GNN的几个问题：

计算复杂度高。原始的Spectral GNN需要对图Laplacian做谱分解，但可以用多项式滤波或类似power iteration的更新方案避免完整特征分解。
Receptive field不同。Spectral GNN原则上可以让图信号传播到整个图，但局部多项式滤波仍然只作用在局部。
可迁移性困难。spectral filter是在某一张图的spectrum上定义的，改变一条边都可能显著改变spectrum，因此把滤波器迁移到新图并不总是容易。