Tech Blog

Pytorch Geometric 系列将结合 pytorch geometric 类库,从原理结合实践来讲解深度图神经网络。和前几期一样,这一系列所有环境已经为大家封装成了 docker image。预制的 image 中既包含了 pytorch 1.8 和对应的 geometric,networkx,Jubyter notebook 还有画图涉及到的 matplotlib 和 hvplot 类库。除此之外,也包含了预下载的 geometric dateset 和所用到的代码。方便大家一键开始 coding。

获取预制 Docker Image

获取docker image 环境的方式,小伙伴们可以通过关注 MyEncyclopedia 公众号并在公众号对话框中回复

docker-geometric 后获得

wechat_docker

关于具体使用我制作的 docker image 这个主题,大家可以再公众号或者同名 B站账号,找到 docker 封装工程依赖 一集, 其中包括如下高级用法

  • 在 VS Code 中运行 Jupyter notebook
  • Remote X Window 主机显示渲染图片
  • docker container 中启动 Jupyter notebook 服务器
  • 共享文件保存文件改动
  • 添加类库后保存镜像改动
bili-docker-adv

如果觉得下载镜像麻烦的小伙伴们,也可以自己准备对应的 python package,再拷贝下方链接的代码直接运行。

https://mydoc.myencyclopedia.top/pub/code/cora-plot

这个系列涉及到的 python package 包括 - torch-geometric

  • networkx

  • hvplot

  • livelossplot

启动 Jupyter Notebook

当预制的 docker image 下载并装载完毕后,通过下面命令来启动 container 内置的 Jupyter notebook 服务器,这样可以使得主机种的浏览器访问到 container 中的服务器。

1
docker run -p 8888:8888 -it env-conda-x jupyter notebook --allow-root --ip 0.0.0.0

内置 Cora 数据集

下一步,我们会通过 geometric 类库加载 cora 数据集。这一步通常来说需要从网上下载,但是预制的 docker image 已经为大家下载好了 planetoid 图数据集。Planetoid 包含了 Cora,Pubmed 和 Citeseer。

因此,加载数据集这一步执行非常快

1
2
3
4
5
6
7
8
9
10
11
12
import numpy as np
import torch
import torch.nn as nn
import torch.nn.functional as F

from torch_geometric.data import Data
from torch_geometric.nn import GATConv
from torch_geometric.datasets import Planetoid
import torch_geometric.transforms as T

name_data = 'Cora'
dataset = Planetoid(root='./data/', name=name_data)

然后我们将 cora 转换成 networkx 格式。networkx 是 python 中一个比较流行的图类库。我们在后面 visualization 中也会利用networkx 的功能。

Cora 有 7 种节点类型,我们将每种节点类型赋予不同颜色,有助于更好 visualization。

1
2
3
4
5
6
7
8
9
10
from torch_geometric.utils import to_networkx
cora = to_networkx(dataset.data)

print(cora.is_directed())

node_classes = dataset.data.y.data.numpy()
print(node_classes)

node_color = ["red","blue","green","yellow","peru","violet","cyan"]
node_label = np.array(list(cora.nodes))

接着,调用 networkx 的 spring_layout 计算每个节点的弹簧布局下的位置,这一步执行会比较耗时。

1
2
3
4
import matplotlib.pyplot as plt
import networkx as nx

pos = nx.layout.spring_layout(cora)

我们首先来看一下 matplotlib 的渲染效果。

1
2
3
4
5
6
7
8
9
10
11
plt.figure(figsize=(16,12))

for i in np.arange(len(np.unique(node_classes))):
node_list = node_label[node_classes == i]
nx.draw_networkx_nodes(cora, pos, nodelist=list(node_list),
node_size=50,
node_color=node_color[i],
alpha=0.8)

nx.draw_networkx_edges(cora, pos,width=1,edge_color="black")
plt.show()

因为 matplotlib 只能画出一张静态图片, 无法做 interaction,也无法动态缩放。因此渲染效果不是特别好,尤其是对于 cora 这种数据量比较大的 graph 尤为显著。

我们看到图片种尽管有七种颜色的节点,但是当中存在的这块密集的点,我们很难看出节点和节点之间的关系。 cora_networkx

我们换一个类库 hvplot,它的渲染和交换效果如下。 hvplot

代码和 matplotlib 大致一致。注意渲染的时候 hvplot 需要将多个图片数据以乘法形式返回,借助 reduce 函数我们将 7 种节点的图相乘,再乘以描绘边的图,呈现出叠加的完整图片。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
import hvplot.networkx as hvnx
options = {
'width': 800,
'height': 1000
}

plt_nodes = []
for i in np.arange(len(np.unique(node_classes))):
nodelist = node_label[node_classes == i]
plt = hvnx.draw_networkx_nodes(cora, pos, nodelist=list(nodelist), node_color=node_color[i], **options)
plt_nodes.append(plt)

plt_edges = hvnx.draw_networkx_edges(cora, pos, arrowstyle='->', edge_width=2, colorbar=True, **options)

import functools
import operator
plt_edges * functools.reduce(operator.mul, plt_nodes)

大家一定都用过神奇的 mathpix 服务,它可以方便大家将数学公式图片转换成 Latex 数学公式。现在 MyEncyclopedia 也推出了手机版和桌面版的 mathpix 服务啦!

使用流程 登录服务

1. 关注 MyEncyclopedia 公众号

2. 在公众号对话中回复 ?

3. 点击 AI在线服务

4. 在非微信浏览器中打开链接

5. 登录

目前提供两种方式。 在国内的小伙伴建议使用微信公众号的方式登录,在国外的朋友们可以使用 Google 账号登录。

5.1 微信登录,按提示在为公众号中输入验证码

5.2 Google登录

使用流程

点击右上方 + 选择图片

等待结果

Output Source 为转换的 Latex 源码

Output Render 为转换源码对应的渲染效果

在学习了一些基本的统计变量生成法之后,这次我们来看看如何生成正态分布。它就是大名鼎鼎的 Box-Muller 方法,Box-Muller 的理解过程可以体会到统计模拟的一些精妙思想。

尝试逆变换方法

我们先尝试通过标准的逆变换方法来生成正态分布。

正态分布的 PDF 表达式为

\[ f_Z(z) = \frac{1}{\sqrt{2 \pi}} \exp\left\{-\frac{z^2}{2}\right\} \]

对应的函数图形是钟形曲线

根据 PDF,其 CDF 的积分形式为 \[ \Phi(x)=\frac{1}{\sqrt{2 \pi}} \int_{-\infty}^{x} e^{-t^{2} / 2} d t \]

和所有 PDF CDF 关系一样,\(\Phi(x)\) 表示 \(f_Z\) 累积到 \(x\) 点的面积。

很不幸的是,\(\Phi(x)\) 无法写出一般数学表达式,因而也无法直接用逆变换方法。

二维映射到一维

我们知道,高维正态分布有特殊的性质:它的每一维的分量都是正态分布;单个维度对于其他维度的条件概率分布也是正态分布。

用图来理解这两条性质就是,对于下图的二维正态分布 $ x = [x_1, x_2]^T $,单独的 \(x_1\)\(x_2\) 都服从一维正态分布。

条件概率 \(p(x_2|x_1 \approx1)\) 的PDF 对应图中的红线,显然也是一维正态分布。

写一段简单的代码验证二维正态分布的单个分量服从正态分布。

代码中,我们用np.random.normal生成了 10000 个服从二维正态分布的 x, y 点,然后我们丢弃 y,只保留 x,并画出 10000 个 x 的分布。

1
2
3
4
5
6
7
8
def plot_normal_1d():
x, _ = np.random.normal(loc=0, scale=1, size=(2, 10000))
import seaborn as sns
sns.distplot(x, hist=True, kde=True, bins=100, color='darkblue',
hist_kws={'edgecolor': 'black'},
kde_kws={'linewidth': 4})
plt.title('PDF Normal 1D from 2D')
plt.show()

Box-Muller 原理

虽然无法直接用逆变换方法生成一维正态分布,但我们却能通过先生成二维的正态分布,利用上面一节的性质,生成一维正态分布。

而 Box-Muller 就是巧妙生成二维正态分布样本点的方法。

首先,我们来看看二维正态分布可以认为是两个维度是独立的,每个维度都是正态分布。此时,其 PDF 可以写成两个一维正态分布 PDF 的乘积。

这种写法表明,二维正态分布仅用一个 r 向量就可以充分表达。注意,r 是向量,不仅有大小还有角度,有两个分量。这两个分量本质上是独立的,这就是 Box-Muller 方法的巧妙之处。也就是,Box-Muller 通过角度和半径大小两个分量的独立性分别单独生成并转换成 (x, y) 对。

角度分量是在 \(2\pi\) 范围均匀采样,这一点比较直觉好理解。

再来看看半径分量 r。我们令 \[ s = {r^2 \over 2} \Longrightarrow r = \sqrt{2s} \]

则 s 服从指数分布 \(\lambda=1\)

不信么?我们不妨来做个模拟实验,下图是模拟 10000次二维正态分布 (x, y) 点后转换成 s 的分布。

模拟和plot 代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
def plot_r_squared():
def gen_normal_samples(n):
x, y = np.random.normal(loc=0, scale=1, size=(2, n))
return x, y

x, y = gen_normal_samples(10000)
s = (x * x + y * y)/2
plot_dist_1d(s, title='PDF $s = {{x^2 + y^2}\over{2}} \sim exp(1)$')


def plot_dist_1d(X, title='PDF '):
import seaborn as sns
plt.rcParams.update({
"text.usetex": True,
"font.family": "sans-serif",
"font.sans-serif": ["Helvetica"]})
sns.distplot(X, hist=True, kde=True, bins=100, color='darkblue',
hist_kws={'edgecolor': 'black'},
kde_kws={'linewidth': 4})
plt.title(title)
plt.show()

确信了 s 符合指数分布,根据指数分布的 PDF,可以推出二维正态 PDF中的 $ e{-r2/2}$ 也符合指数分布,即 \[ s \sim \exp(1) \Longrightarrow e^{-r^2/2} \sim \exp(1) \]

至此,总结一下Box-Muller方法。我们视二维正态分布PDF为独立两部分的乘积,第一部分是在 \(2\pi\) 范围中的均匀分布,代表了二维平面中的角度 \(\theta\),第二部分为 \(\lambda=1\) 的指数分布,代表半径大小。

Box-Muller 方法通过两个服从 [0, 1] 均匀分布的样本 u1和u2,转换成独立的角度和半径样本,具体过程如下

  1. 生成 [0, 1] 的均匀分布 u1,利用逆变换采样方法转换成 exp(1) 样本,此为二维平面点半径 r

  2. 生成 [0, 1] 的均匀分布 u2,乘以 \(2\pi\),即为样本点的角度 \(\theta\)

  3. 将 r 和 \(\theta\) 转换成 x, y 坐标下的点。

理解了整个过程的意义,下面的代码就很直白。

1
2
3
4
5
6
7
8
9
10
def normal_box_muller():
import random
from math import sqrt, log, pi, cos, sin
u1 = random.random()
u2 = random.random()
r = sqrt(-2 * log(u1))
theta = 2 * pi * u2
z0 = r * cos(theta)
z1 = r * sin(theta)
return z0, z1

接下来,我们来看看 Box-Muller 法生成的二维标准正态分布动画吧

拒绝采样极坐标方法

Box-Muller 方法还有一种形式,称为极坐标形式,属于拒绝采样方法。

1. 生成独立的 u, v 和 s

分别生成 [0, 1] 均匀分布 u 和 v。令 \(s = r^2 = u^2 + v^2\)。如果 s = 0或 s ≥ 1,则丢弃 u 和 v ,并尝试另一对 (u , v)。因为 u 和 v 是均匀分布的,并且因为只允许单位圆内的点,所以 s 的值也将均匀分布在开区间 (0, 1) 中。注意,这里的 s 的意义虽然也为半径,但不同于基本方法中的 s。这里 s 取值范围为 (0, 1) ,目的是通过 s 生成指数分布,而基本方法中的 s 取值范围为 [0, +∞],表示二维正态分布 PDF 采样点的半径。复用符号 s 的原因是为了对应维基百科中关于基本方法和极坐标方法的数学描述。

我们用代码来验证 s 服从 (0, 1) 范围上的均匀分布。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
def gen_polar_s():
import random
while True:
u = random.uniform(-1, 1)
v = random.uniform(-1, 1)
s = u * u + v * v
if s >= 1.0 or s == 0.0:
continue
return s


def plot_polar_s():
s = [gen_polar_s() for _ in range(10000) ]
plot_dist_1d(s, title='PDF Polar $s = u^2 + v^2$')

2. 将 u, v, s 转换成 x, y

若将 $s = R^2 uniform(0, 1) $ 看成是基本方法中的 u1,就可以用同样的方式转换成指数分布,用以代表二维PDF的半径。

同时,根据下图,\(\cos \theta\)\(\sin \theta\) 可以直接用 u, v, R 表示出来,并不需要通过三角函数显示计算出 \(\theta\)。有了半径, \(\cos \theta\)\(\sin \theta\) ,则可以直接计算出 x, y 坐标,(下面用 \(z_0, z_1\) 代替 \(x, y\))。

\[ z_{0}=\sqrt{-2 \ln U_{1}} \cos \left(2 \pi U_{2}\right)=\sqrt{-2 \ln s}\left(\frac{u}{\sqrt{s}}\right)=u \cdot \sqrt{\frac{-2 \ln s}{s}} \]

\[ z_{1}=\sqrt{-2 \ln U_{1}} \sin \left(2 \pi U_{2}\right)=\sqrt{-2 \ln s}\left(\frac{v}{\sqrt{s}}\right)=v \cdot \sqrt{\frac{-2 \ln s}{s}} \]

同样,Box-Muller 极坐标方法的代码和公式一致。

1
2
3
4
5
6
7
8
9
10
11
12
def normal_box_muller_polar():
import random
from math import sqrt, log
while True:
u = random.uniform(-1, 1)
v = random.uniform(-1, 1)
s = u * u + v * v
if s >= 1.0 or s == 0.0:
continue
z0 = u * sqrt(-2 * log(s) / s)
z1 = v * sqrt(-2 * log(s) / s)
return z0, z1

拒绝采样的效率

极坐标方法与基本方法的不同之处在于它是一种拒绝采样。因此,它会丢弃一些生成的随机数,但可能比基本方法更快,因为它计算更简单:避免使用昂贵的三角函数,并且在数值上更稳健。极坐标方法丢弃了生成总输入对的 1 − π /4 ≈ 21.46%,即需要 4/ π ≈ 1.2732 个输入随机数,输出一个随机采样。

BiliBili

Lecture 1 - Introduction to Unleashing Novel Data at Scale

  • Why data accessibility is so central to the advancement of knowledge in economics (with some historical background)
  • An overview of the data curation pipeline
    • Step 1: Detect document layouts
    • Step 2: OCR
    • Step 3: Post-processing and database assembly
    • Step 4: Convert information into computable format
  • Why is the material covered in this course useful to social scientists?
    • Why there won’t be an app/commercial product capable of end-to-end processing of social science documents anytime soon
    • Why manual data entry often falls short
    • Why our problems differ from those that are the central focus of computer science and the digital humanities
  • At its core, deep learning is an optimization problem, which economists are well-trained to understand. It would be really unfortunate if we did not take full advantage of the very powerful methods that deep learning offers, which we are well poised to utilize

Lecture 2 - Why Deep Learning?

This post compares rule-based and deep learning-based approaches to data curation. It discusses their requirements and why rule-based approaches often (but not always) produce disappointing results when applied to social science data.

  • An overview of the syllabus (ultimately, the course had a few deviations from the original syllabus, based on student interests; final syllabus posted in the course section of this website)
  • There are two distinct approaches to automated data curation
    • Tell the computer how to process the data by defining a set of rules
    • Let the computer learn how to process the data from empirical examples, using deep learning
  • Overview of rules, how they are used to process image scans and text, why they often fail, and why they sometimes succeed
  • Deep learning, how it contrasts with rule-based approaches, and its requirements
  • Does the noise from rule-based approaches really matter?

Lecture 4 - Convolutional Neural Networks

  • A brief overview of convolutions
  • Benchmark datasets for image classification (following the ConvNent literature requires familiarity with the benchmarks)
  • Image classification with a linear classifier (and its shortcomings)
  • CNN architectures
    • AlexNet
    • VGG
    • GoogLeNet
    • ResNet
    • ResNext

Lecture 5 - Image Classification; Training Neural Nets

This post covers two topics: using CNNs for image classification (a very useful task) and training neural networks in practice. Much of the information about training neural nets is essential to implementing deep learning-based approaches, whether with CNNs or with some other architecture.

Image Classification

  • Loss functions for classification
    • SVM
    • Softmax
  • Deep document classification

Training Neural Nets

  • Activation functions
  • Data pre-processing
  • Initialization
  • Optimization
  • Regularization
  • Batch normalization
  • Dropout
  • Data augmentation
  • Transfer learning
  • Setting hyperparameters
  • Monitoring the learning process

Lecture 6 - Other Computer Vision Problems (Including Object Detection)

This post covers object detection as well as the related problems of semantic segmentation, localization, and instance segmentation. Object detection is core to document image analysis, as it is used to determine the coordinates and classes of different document layout regions. The other problems covered are closely related.

  • Semantic segmentation
  • Localization
  • Object detection
    • Region CNNs
    • Fast R-CNN
    • Faster R-CNN
    • Mask R-CNN
    • Features pyramids
  • Instance segmentation
  • Other frameworks (YOLO)

Lecture 7 - Object Detection in Practice

  • Selecting an object detection model
  • Overview of Detectron2
  • How-to in D2

Lecture 8 - Labeling and Deep Visualization

Labeling

  • Active learning for layout annotation
  • Labeling hacks

Deep visualization

  • Basic visualization approaches
  • Gradient based ascent
  • Deep Dream

Lecture 9 - Generative Adversarial Networks

  • Overview: supervised and unsupervised learning; generative models
  • Generative adversarial networks
  • CycleGAN

Lecture 10 - OCR Architecture

  • Overview of the OCR problem
  • Recurrent neural networks
  • LSTMs
  • Connectionist temporal classification
  • Putting it together

Lecture 11 - OCR and Post-Processing in Practice

This post discusses OCR, both off-the-shelf and how to implement a customized OCR model. It discusses how Layout Parser can be used for end-to-end document image analysis, and provides concrete examples of creating variable domains during post-processing. It also provides an overview of the second half of the knowledge base, which covers NLP.

  • Off-the-shelf OCR

  • Designing customized OCR

  • Putting it altogether (and Layout Parser)

  • Creating variable domains

  • An overview of the second half of the course (NLP)

Lecture 12 - Models of Words

  • Traditional models of words

  • Word2Vec

  • GloVe

  • Evaluation

  • Interpreting word vectors

  • Problems with word vectors

Lecture 13 - Language Modeling and Other Topics in NLP

This post provides an introduction to language modeling, as well as several other important topics: dependency parsing, named entity recognition (NER), and labeling for NLP. Due to time constraints, the course is able to provide only a very brief introduction to topics like dependency parsing and NER, which have traditionally been quite central questions in NLP research.

  • Language Modeling
    • Count based models
    • Bag of words
    • RNN (review)
    • LSTM (review)
  • Dependency parsing
  • Named entity recognition
  • Labeling for NLP

Lecture 14 - Seq2Seq and Machine Translation

Machine translation has pioneered some of the most productive innovations in neural-based NLP and hence is useful to study even for those who care little about machine translation per se. We will focus in particular on seq2seq and attention.

  • Statistical machine translation
  • Neural machine translation

Lecture 15 - Attention is All You Need

This post introduces the Transformer, a seq2seq model based entirely on attention that has transformed NLP. Given the importance of this paper, there are a bunch of very well-done web resources about it, cited in the lecture and below, that I recommend checking out directly (there are others who have much more of a comparative advantage in presenting seminal NLP papers than I do!).

  • A recap of attention

  • The Transformer

    • The encoder
      • Encoder self-attention
      • Positional embeddings
      • Add and normalize
    • The decoder
      • Encoder-decoder attention
      • Decoder self-attention
    • Linear and softmax layers
    • Training

Lecture 16 - Transformer-Based Language Models

This post provides an overview of various Transformer-based language models, discussing their architectures and which are best-suited for different contexts.

  • Overview
  • Contextualized word embeddings
  • Models
    • GPT
    • BERT
    • RoBERTa
    • DistilBERT
    • ALBERT
    • T5
    • GPT2/GPT3
    • Transformers XL
    • XLNet
    • Longformer
    • BigBird
  • Recap and what to use

Lecture 17 - Understanding Transformers, Visualization, and Sentiment Analysis

This post covers a variety of topics around Transformer-based language models: understanding how Transformer attention works, understanding what information is contained in their embeddings, visualizing embeddings, and using Transformer-based models to conduct sentiment analysis.

  • What do Transformer-based models attend to?

  • What’s in an embedding?

  • Visualizing embeddings

  • Sentiment analysis

Lecture 18 - NLP with Noisy Text

  • The Canonical Deep NLP Training Corpus
  • A definition of noise
  • The problem with noise
  • Approaches for denoising

Lecture 19 - Retrieval and Question Answering

  • Reading comprehension
  • Open-domain question answering

Lecture 20 - Zero-Shot and Few-Shot Learning in NLP

  • What it means to learn in just a few shots

  • Zero-shot and few-shot learning in practice

Lecture 21 - Transformers for Computer Vision

  • Transformers for computer vision
  • Transformers for image classification
  • Transformers for object detection

YouTube

BiliBili

在这篇文章中,我们从一道LeetCode 470 题目出发,通过系统地思考,引出拒绝采样(Reject Sampling)的概念,并探索比较三种拒绝采样地解法;接着借助状态转移图来定量计算采样效率;最后,我们利用同样的方法来解一道稍微复杂些的经典抛硬币求期望的统计面试题目。

Leetcode 470 用 Rand7() 实现 Rand10()

已有方法 rand7 可生成 1 到 7 范围内的均匀随机整数,试写一个方法 rand10 生成 1 到 10 范围内的均匀随机整数。

不要使用系统的 Math.random() 方法。

思考

  • rand7()调用次数的 期望值 是多少 ?

  • 你能否尽量少调用 rand7() ?

思维过程

我们已有 rand7() 等概率生成了 [1, 7] 中的数字,我们需要等概率生成 [1, 10] 范围内的数字。第一反应是调用一次rand7() 肯定是不够的,因为覆盖的范围不够。那么,就需要至少2次调用 rand7() 才能生成一次 rand10(),但是还要保证 [1, 10] 的数字生成概率相等,这个是难点。 现在我们先来考虑反问题,给定rand10() 生成 rand7()。这个应该很简单,调用一次 rand10() 得到 [1, 10],如果是 8, 9, 10 ,则丢弃,重新开始,否则返回。想必大家都能想到这个朴素的方法,这种思想就是统计模拟中的拒绝采样(Reject Sampling)。

有了上面反问题的思考,我们可能会想到,rand7() 可以生成 rand5(),覆盖 [1, 5]的范围,如果将区间 [1, 10] 分成两个5个值的区间 [1, 5] 和 [6, 10],那么 rand7() 可以通过先等概率选择区间 [1, 5] 或 [6, 10],再通过rand7() 生成 rand5()就可以了。这个问题就等价于先用 rand7() 生成 rand2(),决定了 [1, 5] 还是 [6, 10],再通过rand7() 生成 rand5() 。

解法一:rand2() + rand5()

我们来实现这种解法。下图为调用两次 rand7() 生成 rand10 数值的映射关系:横轴表示第一次调用,1,2,3决定选择区间 [1, 5] ,4,5,6选择区间 [6, 10]。灰色部分表示结果丢弃,重新开始 (注,若第一次得到7无需再次调用 rand7())。

有了上图,我们很容易写出如下 AC 代码。

{linenos
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# AC
# Runtime: 408 ms, faster than 23.80% of Python3 online submissions for Implement Rand10() Using Rand7().
# Memory Usage: 16.7 MB, less than 90.76% of Python3 online submissions for Implement Rand10() Using Rand7().
class Solution:
def rand10(self):
while True:
a = rand7()
if a <= 3:
b = rand7()
if b <= 5:
return b
elif a <= 6:
b = rand7()
if b <= 5:
return b + 5

标准解法:rand49()

从提交的结果来看,第一种解法慢于多数解法。原因是我们的调用 rand7() 的采样效率比较低,第一次有 1/7 的概率结果丢弃,第二次有 2/7的概率被丢弃。

如何在第一种解法的基础上提高采样效率呢?直觉告诉我们一种做法是降低上述 7x7 表格中灰色格子的面积。此时,会想到我们通过两次 rand7() 已经构建出来 rand49()了,那么再生成 rand10() 也规约成基本问题了。

下图为 rand49() 和 rand10() 的数字对应关系。

实现代码比较简单。注意,while True 可以去掉,用递归来代替。

{linenos
1
2
3
4
5
6
7
8
9
10
11
# AC
# Runtime: 376 ms, faster than 54.71% of Python3 online submissions for Implement Rand10() Using Rand7().
# Memory Usage: 16.9 MB, less than 38.54% of Python3 online submissions for Implement Rand10() Using Rand7().
class Solution:
def rand10(self):
while True:
a, b = rand7(), rand7()
num = (a - 1) * 7 + b
if num <= 40:
return num % 10 + 1

更快的做法

上面的提交结果发现标准解法在运行时间上有了不少提高,处于中等位置。我们继续思考,看看能否再提高采样效率。

观察发现,rand49() 有 9/49 的概率,生成的值被丢弃,原因是 [41, 49] 只有 9 个数,不足10个。倘若此时能够将这种状态保持下去,那么只需再调用一次 rand7() 而不是从新开始情况下至少调用两次 rand7(), 就可以得到 rand10()了。也就是说,当 rand49() 生成了 [41, 49] 范围内的数的话等价于我们先调用了一次 rand9(),那么依样画葫芦,我们接着调用 rand7() 得到了 rand63()。63 分成了6个10个值的区间后,剩余 3 个数。此时,又等价于 rand3(),循环往复,调用了 rand7() 得到了 rand21(),最后若rand21() 不幸得到21,等价于 rand1(),此时似乎我们走投无路,只能回到最初的状态,一切从头再来了。

改进算法代码如下。注意这次击败了 92.7%的提交。

{linenos
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# AC
# Runtime: 344 ms, faster than 92.72% of Python3 online submissions for Implement Rand10() Using Rand7().
# Memory Usage: 16.7 MB, less than 90.76% of Python3 online submissions for Implement Rand10() Using Rand7().
class Solution:
def rand10(self):
while True:
a, b = rand7(), rand7()
num = (a - 1) * 7 + b
if num <= 40: return num % 10 + 1
a = num - 40
b = rand7()
num = (a - 1) * 7 + b
if num <= 60: return num % 10 + 1
a = num - 60
b = rand7()
num = (a - 1) * 7 + b
if num <= 20: return num % 10 + 1

采样效率计算

通过代码提交的结果和大致的分析,我们已经知道三个解法在采样效率依次变得更快。现在我们来定量计算这三个解法。

我们考虑生成一个 rand10() 平均需要调用多少次 rand7(),作为采样效率的标准。

一种思路是可以通过模拟方法,即将上述每个解法模拟多次,然后用总的 rand7() 调用次数除以 rand10() 的生成次数即可。下面以解法三为例写出代码

{linenos
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
# The rand7() API is already defined for you.
rand7_c = 0
rand10_c = 0

def rand7():
global rand7_c
rand7_c += 1
import random
return random.randint(1, 7)

def rand10():
global rand10_c
rand10_c += 1
while True:
a, b = rand7(), rand7()
num = (a - 1) * 7 + b
if num <= 40: return num % 10 + 1
a = num - 40 # [1, 9]
b = rand7()
num = (a - 1) * 7 + b # [1, 63]
if num <= 60: return num % 10 + 1
a = num - 60 # [1, 3]
b = rand7()
num = (a - 1) * 7 + b # [1, 21]
if num <= 20: return num % 10 + 1

if __name__ == '__main__':
while True:
rand10()
print(f'{rand10_c} {round(rand7_c/rand10_c, 2)}')

运行代码,发现解法三的采样效率稳定在 2.19。

采样效率精确计算

计算解法二

为了精确计算三个解法的采样效率,我们通过代码得到对应的状态转移图来帮助计算。

例如,解法一可以对应到下图:初始状态 Start 节点中的 +2 表示经过此节点会产生 2次 rand7() 的代价。从 Start 节点有 40/49 的概率到达被接受状态 AC,有 9/49 概率到达拒绝状态 REJ。REJ 需要从头开始,则用虚线表示重新回到 Start节点,也就是说 REJ 的代价等价于 Start。注意,从某个节点出发的所有边之和为1。

有了上述状态转移关系图,我们令初始状态的平均代价为 \(C_2\),则可以写成递归表达式,因为其中 REJ 的代价就是 \(C_2\),即

\[ C_2 = 2 + (\frac{40}{49}\cdot0 + \frac{9}{49} C_2) \]

解得 \(C_2\)

\[ C_2 = 2.45 \]

计算解法一

同样的,对于另外两种解法,虽然略微复杂,也可以用同样的方法求得。

解法一的状态转移图为

递归方程表达式为

\[ C_1 = 1+\frac{3}{7} \cdot (1+\frac{5}{7} \cdot 0 + \frac{2}{7} \cdot C_1) \cdot2+ \frac{1}{7} \cdot (C_1) \]

解得 \(C_1\)

\[ C_1 = \frac{91}{30} \approx 3.03 \]

计算解法三

最快的解法三状态转移图为

递归方程表达式为

\[ C_3 = 2+\frac{40}{49} \cdot 0 + \frac{9}{49} (1+\frac{60}{63} \cdot 0 + \frac{3}{63} \cdot (1+\frac{20}{21} \cdot 0 + \frac{1}{21} \cdot C_3)) \]

解得 \(C_3\) \[ C_3 = \frac{329}{150} \approx 2.193 \]

至此总结一下,三个解法的平均代价为 \[ C_1 \approx 3.03 > C_2 = 2.45 > C_3 \approx 2.193 \] 这些值和我们通过模拟得到的结果一致。

稍难些的经典概率求期望题目

至此,LeetCode 470 我们已经分析透彻。现在,我们已经可以很熟练的将此类拒绝采样的问题转变成有概率的状态转移图,再写成递推公式去求平均采样的代价(即期望)。这里,如果大家感兴趣的话不妨再来看一道略微深入的经典统计概率求期望的题目。

问题:给定一枚抛正反面概率一样的硬币,求连续抛硬币直到两个正面(正面记为H,两个正面HH)的平均次数。例如:HTTHH是一个连续次数为5的第一次出现HH的序列。

分析问题画出状态转移图:我们令初始状态下得到第一个HH的平均长度记为 S,那么下一次抛硬币有 1/2 机会是 T,此时状态等价于初始状态,另有 1/2 机会是 H,我们记这个状态下第一次遇见HH的平均长度为 H(下图蓝色节点)。从此蓝色节点出发,当下一枚硬币是H则结束,是T是返回初始状态。于是构建出下图。

这个问题稍微复杂的地方在于我们有两个未知状态互相依赖,但问题的本质和方法是一样的,分别从 S 和 H 出发考虑状态的概率转移,可以写成如下两个方程式:

\[ \left\{ \begin{array}{c} S =&\frac{1}{2} \cdot(1+H) + \frac{1}{2} \cdot(1+S) \\ H =&\frac{1}{2} \cdot 1 + \frac{1}{2} \cdot(1+S) \end{array} \right. \]

解得

\[ \left\{ \begin{array}{c} H= 4 \\ S = 6 \end{array} \right. \]

因此,平均下来,需要6次抛硬币才能得到 HH,这个是否和你直觉的猜测一致呢?

这个问题还可以有另外一问,可以作为思考题让大家来练习一下:第一次得到 HT 的平均次数是多少?这个是否和 HH 一样呢?

本篇是深度强化学习动手系列文章,自MyEncyclopedia公众号文章深度强化学习之:DQN训练超级玛丽闯关发布后收到不少关注和反馈,这一期,让我们实现目前主流深度强化学习算法PPO来打另一个红白机经典游戏1942。

相关文章链接如下:

强化学习开源环境集

视频论文解读:PPO算法

视频论文解读:组合优化的强化学习方法

解读TRPO论文,深度强化学习结合传统优化方法

解读深度强化学习基石论文:函数近似的策略梯度方法

NES 1942 环境安装

红白机游戏环境可以由OpenAI Retro来模拟,OpenAI Retro还在 Gym 集成了其他的经典游戏环境,包括Atari 2600,GBA,SNES等。

不过,受到版权原因,除了一些基本的rom,大部分游戏需要自行获取rom。

环境准备部分相关代码如下

1
pip install gym-retro
1
python -m retro.import /path/to/your/ROMs/directory/

OpenAI Gym 输入动作类型

在创建 retro 环境时,可以在retro.make中通过参数use_restricted_actions指定 action space,即按键的配置。

1
env = retro.make(game='1942-Nes', use_restricted_actions=retro.Actions.FILTERED)

可选参数如下,FILTERED,DISCRETE和MULTI_DISCRETE 都可以指定过滤的动作,过滤动作需要通过配置文件加载。

1
2
3
4
5
6
7
8
class Actions(Enum):
"""
Different settings for the action space of the environment
"""
ALL = 0 #: MultiBinary action space with no filtered actions
FILTERED = 1 #: MultiBinary action space with invalid or not allowed actions filtered out
DISCRETE = 2 #: Discrete action space for filtered actions
MULTI_DISCRETE = 3 #: MultiDiscete action space for filtered actions

DISCRETE和MULTI_DISCRETE 是 Gym 里的 Action概念,它们的基类都是gym.spaces.Space,可以通过 sample()方法采样,下面具体一一介绍。

  • Discrete:对应一维离散空间,例如,Discrete(n=4) 表示 [0, 3] 范围的整数。
1
2
3
from gym.spaces import Discrete
space = Discrete(4)
print(space.sample())

输出是

1
3
  • Box:对应多维连续空间,每一维的范围可以用 [low,high] 指定。 举例,Box(low=-1.0, high=2, shape=(3, 4,), dtype=np.float32) 表示 shape 是 [3, 4],每个范围在 [-1, 2] 的float32型 tensor。
1
2
3
4
from gym.spaces import Box
import numpy as np
space = Box(low=-1.0, high=2.0, shape=(3, 4), dtype=np.float32)
print(space.sample())

输出是

1
2
3
[[-0.7538084   0.96901214  0.38641307 -0.05045208]
[-0.85486996 1.3516271 0.3222616 1.2540635 ]
[-0.29908678 -0.8970335 1.4869047 0.7007356 ]]

  • MultiBinary: 0或1的多维离散空间。例如,MultiBinary([3,2]) 表示 shape 是3x2的0或1的tensor。
    1
    2
    3
    from gym.spaces import MultiBinary
    space = MultiBinary([3,2])
    print(space.sample())

输出是

1
2
3
[[1 0]
[1 1]
[0 0]]
  • MultiDiscrete:多维整型离散空间。例如,MultiDiscrete([5,2,2]) 表示三维Discrete空间,第一维范围在 [0-4],第二,三维范围在[0-1]。
1
2
3
from gym.spaces import MultiDiscrete
space = MultiDiscrete([5,2,2])
print(space.sample())

输出是

1
[2 1 0]
  • Tuple:组合成 tuple 复合空间。举例来说,可以将 Box,Discrete,Discrete组成tuple 空间:Tuple(spaces=(Box(low=-1.0, high=1.0, shape=(3,), dtype=np.float32), Discrete(n=3), Discrete(n=2)))
1
2
3
4
from gym.spaces import *
import numpy as np
space = Tuple(spaces=(Box(low=-1.0, high=1.0, shape=(3,), dtype=np.float32), Discrete(n=3), Discrete(n=2)))
print(space.sample())

输出是

1
2
(array([ 0.22640526,  0.75286865, -0.6309239 ], dtype=float32), 0, 1)

  • Dict:组合成有名字的复合空间。例如,Dict({'position':Discrete(2), 'velocity':Discrete(3)})
    1
    2
    3
    from gym.spaces import *
    space = Dict({'position':Discrete(2), 'velocity':Discrete(3)})
    print(space.sample())

输出是

1
OrderedDict([('position', 1), ('velocity', 1)])

NES 1942 动作空间配置

了解了 gym/retro 的动作空间,我们来看看1942的默认动作空间

1
2
env = retro.make(game='1942-Nes')
print("The size of action is: ", env.action_space.shape)

1
The size of action is:  (9,)

表示有9个 Discrete 动作,包括 start, select这些控制键。

从训练1942角度来说,我们希望指定最少的有效动作取得最好的成绩。根据经验,我们知道这个游戏最重要的键是4个方向加上 fire 键。限定游戏动作空间,官方的做法是在创建游戏环境时,指定预先生成的动作输入配置文件。但是这个方式相对麻烦,我们采用了直接指定按键的二进制表示来达到同样的目的,此时,需要设置 use_restricted_actions=retro.Actions.FILTERED。

下面的代码限制了6种按键,并随机play。

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
action_list = [
# No Operation
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
# Left
[0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
# Right
[0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
# Down
[0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
# Up
[0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
# B
[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
]

def random_play(env, action_list, sleep_seconds=0.01):
env.viewer = None
state = env.reset()
score = 0
for j in range(10000):
env.render()
time.sleep(sleep_seconds)
action = np.random.randint(len(action_list))

next_state, reward, done, _ = env.step(action_list[action])
state = next_state
score += reward
if done:
print("Episode Score: ", score)
env.reset()
break

env = retro.make(game='1942-Nes', use_restricted_actions=retro.Actions.FILTERED)
random_play(env, action_list)

来看看其游戏效果,全随机死的还是比较快。

图像输入处理

一般对于通过屏幕像素作为输入的RL end-to-end训练来说,对图像做预处理很关键。因为原始图像较大,一方面我们希望能尽量压缩图像到比较小的tensor,另一方面又要保证关键信息不丢失,比如子弹的图像不能因为图片缩小而消失。另外的一个通用技巧是将多个连续的frame合并起来组成立体的frame,这样可以有效表示连贯动作。

下面的代码通过 pipeline 将游戏每帧原始图像从shape (224, 240, 3) 转换成 (4, 84, 84),也就是原始的 width=224,height=240,rgb=3转换成 width=84,height=240,stack_size=4的黑白图像。具体 pipeline为

  1. MaxAndSkipEnv:每两帧过滤一帧图像,减少数据量。

  2. FrameDownSample:down sample 图像到指定小分辨率 84x84,并从彩色降到黑白。

  3. FrameBuffer:合并连续的4帧,形成 (4, 84, 84) 的图像输入

1
2
3
4
5
6
7
def build_env():
env = retro.make(game='1942-Nes', use_restricted_actions=retro.Actions.FILTERED)
env = MaxAndSkipEnv(env, skip=2)
env = FrameDownSample(env, (1, -1, -1, 1))
env = FrameBuffer(env, 4)
env.seed(0)
return env

观察图像维度变换

1
2
3
4
5
env = retro.make(game='1942-Nes', use_restricted_actions=retro.Actions.FILTERED)
print("Initial shape: ", env.observation_space.shape)

env = build_env(env)
print("Processed shape: ", env.observation_space.shape)

确保shape 从 (224, 240, 3) 转换成 (4, 84, 84)

1
2
Initial shape:  (224, 240, 3)
Processed shape: (4, 84, 84)

FrameDownSample实现如下,我们使用了 cv2 类库来完成黑白化和图像缩放

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
class FrameDownSample(ObservationWrapper):
def __init__(self, env, exclude, width=84, height=84):
super(FrameDownSample, self).__init__(env)
self.exclude = exclude
self.observation_space = Box(low=0,
high=255,
shape=(width, height, 1),
dtype=np.uint8)
self._width = width
self._height = height

def observation(self, observation):
# convert image to gray scale
screen = cv2.cvtColor(observation, cv2.COLOR_RGB2GRAY)

# crop screen [up: down, left: right]
screen = screen[self.exclude[0]:self.exclude[2], self.exclude[3]:self.exclude[1]]

# to float, and normalized
screen = np.ascontiguousarray(screen, dtype=np.float32) / 255

# resize image
screen = cv2.resize(screen, (self._width, self._height), interpolation=cv2.INTER_AREA)
return screen

MaxAndSkipEnv,每两帧过滤一帧

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
class MaxAndSkipEnv(Wrapper):
def __init__(self, env=None, skip=4):
super(MaxAndSkipEnv, self).__init__(env)
self._obs_buffer = deque(maxlen=2)
self._skip = skip

def step(self, action):
total_reward = 0.0
done = None
for _ in range(self._skip):
obs, reward, done, info = self.env.step(action)
self._obs_buffer.append(obs)
total_reward += reward
if done:
break
max_frame = np.max(np.stack(self._obs_buffer), axis=0)
return max_frame, total_reward, done, info

def reset(self):
self._obs_buffer.clear()
obs = self.env.reset()
self._obs_buffer.append(obs)
return obs

FrameBuffer,将最近的4帧合并起来

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
class FrameBuffer(ObservationWrapper):
def __init__(self, env, num_steps, dtype=np.float32):
super(FrameBuffer, self).__init__(env)
obs_space = env.observation_space
self._dtype = dtype
self.observation_space = Box(low=0, high=255, shape=(num_steps, obs_space.shape[0], obs_space.shape[1]), dtype=self._dtype)

def reset(self):
frame = self.env.reset()
self.buffer = np.stack(arrays=[frame, frame, frame, frame])
return self.buffer

def observation(self, observation):
self.buffer[:-1] = self.buffer[1:]
self.buffer[-1] = observation
return self.buffer

最后,visualize 处理后的图像,同样还是在随机play中,确保关键信息不丢失

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
def random_play_preprocessed(env, action_list, sleep_seconds=0.01):
import matplotlib.pyplot as plt

env.viewer = None
state = env.reset()
score = 0
for j in range(10000):
time.sleep(sleep_seconds)
action = np.random.randint(len(action_list))

plt.imshow(state[-1], cmap="gray")
plt.title('Pre Processed image')
plt.pause(sleep_seconds)

next_state, reward, done, _ = env.step(action_list[action])
state = next_state
score += reward
if done:
print("Episode Score: ", score)
env.reset()
break

matplotlib 动画输出

CNN Actor & Critic

Actor 和 Critic 模型相同,输入是 (4, 84, 84) 的图像,输出是 [0, 5] 的action index。

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
class Actor(nn.Module):
def __init__(self, input_shape, num_actions):
super(Actor, self).__init__()
self.input_shape = input_shape
self.num_actions = num_actions

self.features = nn.Sequential(
nn.Conv2d(input_shape[0], 32, kernel_size=8, stride=4),
nn.ReLU(),
nn.Conv2d(32, 64, kernel_size=4, stride=2),
nn.ReLU(),
nn.Conv2d(64, 64, kernel_size=3, stride=1),
nn.ReLU()
)

self.fc = nn.Sequential(
nn.Linear(self.feature_size(), 512),
nn.ReLU(),
nn.Linear(512, self.num_actions),
nn.Softmax(dim=1)
)

def forward(self, x):
x = self.features(x)
x = x.view(x.size(0), -1)
x = self.fc(x)
dist = Categorical(x)
return dist

PPO核心代码

先计算 \(r_t(\theta)\),这里采用了一个技巧,对 \(\pi_\theta\) 取 log,相减再取 exp,这样可以增强数值稳定性。

1
2
3
4
dist = self.actor_net(state)
new_log_probs = dist.log_prob(action)
ratio = (new_log_probs - old_log_probs).exp()
surr1 = ratio * advantage

surr1 对应PPO论文中的 \(L^{CPI}\)

然后计算 surr2,对应 \(L^{CLIP}\) 中的 clip 部分,clip可以由 torch.clamp 函数实现。\(L^{CLIP}\) 则对应 actor_loss。

1
2
surr2 = torch.clamp(ratio, 1.0 - self.clip_param, 1.0 + self.clip_param) * advantage
actor_loss = - torch.min(surr1, surr2).mean()

最后,计算总的 loss \(L_t^{CLIP+VF+S}\),包括 actor_loss,critic_loss 和 policy的 entropy。

1
2
3
4
entropy = dist.entropy().mean()

critic_loss = (return_ - value).pow(2).mean()
loss = actor_loss + 0.5 * critic_loss - 0.001 * entropy

上述完整代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
for _ in range(self.ppo_epoch):
for state, action, old_log_probs, return_, advantage in sample_batch():
dist = self.actor_net(state)
value = self.critic_net(state)

entropy = dist.entropy().mean()
new_log_probs = dist.log_prob(action)

ratio = (new_log_probs - old_log_probs).exp()
surr1 = ratio * advantage
surr2 = torch.clamp(ratio, 1.0 - self.clip_param, 1.0 + self.clip_param) * advantage

actor_loss = - torch.min(surr1, surr2).mean()
critic_loss = (return_ - value).pow(2).mean()

loss = actor_loss + 0.5 * critic_loss - 0.001 * entropy

# Minimize the loss
self.actor_optimizer.zero_grad()
self.critic_optimizer.zero_grad()
loss.backward()
self.actor_optimizer.step()
self.critic_optimizer.step()

补充一下 GAE 的计算,advantage 根据公式

可以转换成如下代码

1
2
3
4
5
6
7
8
9
def compute_gae(self, next_value):
gae = 0
returns = []
values = self.values + [next_value]
for step in reversed(range(len(self.rewards))):
delta = self.rewards[step] + self.gamma * values[step + 1] * self.masks[step] - values[step]
gae = delta + self.gamma * self.tau * self.masks[step] * gae
returns.insert(0, gae + values[step])
return returns

外层 Training 代码

外层调用代码基于随机 play 的逻辑,agent.act()封装了采样和 forward prop,agent.step() 则封装了 backprop 和参数学习迭代的逻辑。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
for i_episode in range(start_epoch + 1, n_episodes + 1):
state = env.reset()
score = 0
timestamp = 0

while timestamp < 10000:
action, log_prob, value = agent.act(state)
next_state, reward, done, info = env.step(action_list[action])
score += reward
timestamp += 1

agent.step(state, action, value, log_prob, reward, done, next_state)
if done:
break
else:
state = next_state

训练结果

让我们来看看学习的效果吧,注意我们的飞机学到了一些关键的技巧,躲避子弹;飞到角落尽快击毙敌机;一定程度预测敌机出现的位置并预先走到位置。

导读:极大似然估计(MLE) 是统计机器学习中最基本的概念,但是能真正全面深入地理解它的性质和背后和其他基本理论的关系不是件容易的事情。极大似然估计和以下概念都有着紧密的联系:随机变量,无偏性质(unbiasedness),一致估计(consistent),asymptotic normality,最优化(optimization),Fisher Information,MAP(最大后验估计),KL-Divergence,sufficient statistics等。在众多阐述 MLE 的文章或者课程中,总体来说都比较抽象,注重公式推导。本系列文章受3blue1brown 可视化教学的启发,坚持从第一性原理出发,通过数学原理结合模拟和动画,深入浅出地让读者理解极大似然估计。

相关链接

抛硬币问题

我们来思考这个老套问题,考虑手上有一枚硬币,旋转(抛)硬币得到正反面的概率固定(令正面概率为\(\theta^{\star}\))但未知,我们如何能通过实验推测出 \(\theta^{\star}\)

朴素的想法是,不断尝试抛硬币,随着次数 n 的增多,正面的比例会趋近于 \(\theta^{\star}\)

对应到数学形式上,令我们对于 \(\theta^{\star}\) 的估计为 \(\hat{\theta}_{n}\),则希望 \[ \hat{\theta}_n = {n_{head} \over n} \to \theta^{\star} \text{ as n } \to \infty \]

模拟试验代码

假设我们尝试了n次,每次的结果为 \(x_i\)\(x_i\)为1(正面) 或 0(反面)。比如试了三次的结果是 [1, 0, 1],则 \(x_1=1, x_2=0, x_3=1\)。一般,我们将观察到的数据写成向量形式

\[X=[x_1, x_2, x_3]^T=[1, 0, 1]^{T}\]

我们知道硬币的正反结果符合伯努利分布,也就是 \[ \begin{align*} P_{ber}(x;\theta) = \left\lbrace \begin{array}{r@{}l} \theta &\text{ if x=1} \\ 1-\theta &\text{ if x=0} \end{array} \right. \end{align*} \]

因为 x 只有0,1两种取值,因此上式也可以写成等价如下的不含条件分支的形式 \[ P_{ber} = \theta^x \cdot (1-\theta)^x \]

假设 \(\theta^{\star} = 0.7\),如果做 n=10 次试验,结果应该比较接近7个1,3个0。

下面我们来模拟一下 n=10,看看结果如何。

下面代码的实现上我们直接使用了pytorch 内置的 bernoulli 函数生成 n 个随机变量实例

1
2
3
4
5
def gen_coins(theta, n=1):
import torch
theta_vec = torch.tensor(n*[theta])
random_values = torch.bernoulli(theta_vec)
return random_values

让我们来做三次 n=10 的试验

1
2
3
4
5
6
for i in range(3):
coins = gen_coins(theta=0.7, n=10)
print(f'trial {i}')
print(f'head #: {sum(coins)}')
print(f'tail #: {sum(1-coins)}')
print()

能发现 7个1,3个0 确实是比较可能的结果。

1
2
3
4
5
6
7
8
9
10
11
trial 0
head #: 7.0
tail #: 3.0

trial 1
head #: 9.0
tail #: 1.0

trial 2
head #: 7.0
tail #: 3.0

生成概率

直觉告诉我们,当 \(\theta^{\star}=0.7\) 时,根据 $P_{ber}(x;) $,7个1,3个0 出现的概率应该是最大,6个1,4个0 或者 8个1,2个0 这两种情况出现概率稍小,其他的情况概率更小。通过基本概率和伯努利公式,重复 n 次试验 1和0出现的概率可以由下面公式算出。(注:7个1,3个0不是单一事件,需要乘以组合数算出实际概率)

\[ P_{X} = \theta^{heads} \cdot (1-\theta)^{tails} \cdot {n \choose heads} \]

P(X)
head=0 0.000006
head=1 0.000138
head=2 0.000032
head=3 0.001447
head=4 0.036757
head=5 0.102919
head=6 0.200121
head=7 0.266828
head=8 0.233474
head=9 0.121061
head=10 0.028248

画出图看的很明显,1出现7次的概率确实最大。

回到我们的问题,我们先假定 \(\theta^{\star} = 0.7\) 的硬币做 n=10 次试验的结果就是 7个1,3个0,或者具体序列为 [1, 0, 0, 1, 0, 1, 1, 1, 1, 1]。那么我们希望按照某种方法推测的估计值 \(\hat\theta\) 也为 0.7。

若将这个方法也记做 \(\hat\theta\),它是\(X\) 的函数,即 \(\hat\theta(X=[1, 0, 0, 1, 0, 1, 1, 1, 1, 1]^T)=0.7\)

我们如何构建这个方法呢?很显然,\(X\) 中 1 的个数就可以胜任,\(\hat\theta=\bar X\)。这个方式确实是正确的,后面的文章我们也会证明它是MLE在伯努利分布参数估计时的计算方法。

但是伯努利分布参数估计的问题中是最简单的情况,背后对应的更一般的问题是:假设我们知道某个过程或者实验生成了某种分布 P,但是不知道它的参数 \(\theta\),如何能通过反复的试验来推断 \(\theta\),同时,我们希望随着试验次数的增多,\(\hat\theta\) 能逼近 \(\theta\)

由于过程是有随机性,试验结果 \(X\) 并不能确定一定是从 \(\hat\theta\) 生成的,因此我们需要对所有 \(\theta\) 打分。对于抛硬币试验来说,我们穷举所有在 [0, 1] 范围内的 \(\theta\),定义它的打分函数 \(f(\theta)\),并且希望我们定义的 \(f(\theta;X=[1, 0, 0, 1, 0, 1, 1, 1, 1, 1]^T)\)\(\theta=0.7\) 时得分最高。推广到一般场景,有如下性质 \[ f(\theta^\star;X) >= f(\theta;X) \]

如此,我们将推测参数问题转换成了优化问题 \[ \hat\theta = \theta^{\star} = \operatorname{argmax}_{\theta} f(\theta; X) = 0.7 \]

朴素方法

一种朴素的想法是,由于 \(\theta^\star=0.7\),因此我们每次的结果应该稍微偏向 1,如果出现了 1,就记0.7分,出现了0,记0.3分,那么我们可以用10个结果的总分来定义总得分,即最大化函数

\[ \begin{equation*} \begin{aligned} &\operatorname{argmax}_{\theta} f(\theta) \\ =& \operatorname{argmax}_{\theta} P(x_1) + P(x_2) + ... + P(x_n) \\ =& \operatorname{argmax}_{\theta} P(x_1|\theta) + P(x_2|\theta) + ... + P(x_n|\theta) \\ =& \operatorname{argmax}_{\theta} \sum P(x_i|\theta) \\ \end{aligned} \end{equation*} \]

很可惜,我们定义的 f 并不符合 \(\theta=0.7\) 时取到最大的原则。下面画出了 \(\theta\) 在 [0, 1] 范围内 f 值,X 固定为 [1, 0, 0, 1, 0, 1, 1, 1, 1, 1]。显然,极值在 0.5 左右。

这种对于观察到的变量实例在整个参数空间打分的方法是最大似然方法的雏形。我们将每次试验结果对于不同 \(\theta\) 的打分就是似然函数的概念。

伯努利似然函数(Likelihood)

伯努利单个结果的似然函数 \(l(\theta)\) 视为 \(\theta\) 的函数,x视为给定值,它等价于概率质量函数 PMF

\[ l(\theta|x) = \theta^x \cdot (1-\theta)^x \]

极大似然估计(MLE)

有了单个结果的似然函数,我们如何定义 \(f(\theta)\) 呢?我们定义的 \(f(\theta)\) 需要满足,在 \(\theta^\star=0.7\)\(n=10\) 的情况下,试验最有可能的结果是 7 个1,3个0,此时 f 需要在 \(\theta=0.7\) 时取到最大值。

极大似然估计(MLE) 为我们定义了合理的 \(f(\theta)\) ,和朴素的想法类似,但是这次用单个结果的似然函数连乘而非连加 \[ L(\theta|X) = l(\theta|x_1) \cdot l(\theta|x_2) \cdot ...l(\theta|x_n) = \prod l(\theta|x_i) \]

我们再来看一下当 $X=[1, 0, 0, 1, 0, 1, 1, 1, 1, 1] $ 时 \(L\)\(\theta\) 空间的取值情况,果然,MLE 能在 0.7时取到最大值。

对数似然函数

最大似然函数 $_{} L() $ 能让我们找到最可能的 \(\theta\),但现实中,我们一般采用最大其 log 的形式。

\[ \begin{equation*} \begin{aligned} &\operatorname{argmax}_{\theta} \log L(\theta|X) \\ =& \operatorname{argmax}_{\theta} \log [l(\theta|x_1) \cdot l(\theta|x_2) \cdot ... \cdot l(\theta|x_n)] \\ =& \operatorname{argmax}_{\theta} \log l(\theta|x_1) + \log l(\theta|x_2) \cdot ... + \log l(\theta|x_n) \end{aligned} \end{equation*} \]

理论能证明,最大对数似然函数得到的极值等价于最大似然函数。但这么做有什么额外好处呢?

我们先将对数似然函数画出来

它的极大值也在 0.7,但是我们发现对数似然函数是个 concave 函数。在优化领域,最大化 concave 函数或者最小化 convex 函数可以有非常高效的解法。再仔细看之前的似然函数,它并不是一个 concave 函数。另一个非常重要的好处是,随着 n 的增大,连乘会导致浮点数 underflow,而单个点的对数似然函数的和的形式就不会有这个问题。

Pytorch MLE 代码

就让我们来实践一下,通过 pytorch 梯度上升来找到极值点。

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
from stats.coin import gen_coins
from collections import deque


def train(num_head: int, num_tail: int) -> float:
import torch
theta = torch.tensor(0.5, requires_grad=True)

recent = deque(3*[100], maxlen=3)

lr = 0.00001
for iter in range(2000):
loss = -(num_head * torch.log(theta) + num_tail * torch.log(1 - theta))
loss.backward()
with torch.no_grad():
theta -= lr * theta.grad
# print(f'{iter}: {theta}, {theta.grad}')
recent.append(theta.grad.item())
if all(map(lambda x: abs(x) < 1, recent)):
break
theta.grad.zero_()
return theta.item()


if __name__ == '__main__':
data = gen_coins(0.6, n=200)

num_head = (data.detach() == 1).sum().item()
num_tail = (data.detach() == 0).sum().item()

print(num_head, num_tail)
print(train(num_head, num_tail))

一点需要说明的是,在迭代过程中,我们保存最后三个导数的值,当最新的三个导数都很小时就退出迭代。

1
if all(map(lambda x: abs(x) < 1, recent))

运行代码,能发现最大化对数似然函数能很稳定的找到 \(\theta\)

现在大家对于伯努利MLE有了一定了解,接着,我们来思考一下最大化似然函数方法是否随着观察次数的增多能不断逼近真实的 \(\theta^\star\)呢?

MLE \(\theta\) 估计的收敛性

\(\theta^\star=0.7\) 的情况下,我们来这样做试验,第一次做 n=1生成观察数据 \(X_{1}\),第二次做 n=2生成观察数据 \(X_{2}\) \[ X_1,X_2, X_3, ..., X_N \] 对于每个数据集 \(X_i\) 通过最大似然方法求得估计的 \(\hat\theta\) \[ \hat\theta_1=MLE(X_1), \hat\theta_2=MLE(X_2), ..., \hat\theta_N=MLE(X_N) \] 将这些 \(\hat\theta_i\) 画出来,可以看到,随着 \(n \to \infty\)\(\hat\theta_i \to \theta^\star=0.7\)

换一个角度来看一下,我们将 \(\hat\theta_i\) 数列按照顺序,离散化后再归一化比例,如下图画出来,红色的柱代表了最新的值 \(\hat\theta\)。可以发现,初始时候,\(\hat\theta\) 在较远离 0.7 的地方出现,随着 n 的增大,出现的位置比较接近 0.7。

MLE \(\theta\) 估计的偏差和方差

我们已经知道 MLE 方法可以通过观察数据推测出最有可能的 \(\hat\theta\),由于观察数据 \(X\) 是伯努利过程产生的,具有随机性,那么 \(\hat\theta\) 可以看成是 \(\theta^\star\) 的随机变量。我们已经通过上面的试验知道随着试验次数的增大,我们的估计会越来越逼近真实值,现在的问题是对于固定的n\(\hat\theta\) 的方差是多少,它的均值是否是无偏的呢?

带着这样的疑问,我们现在做如下试验:

固定 n=10,重复做实验,画出随着次数增多 \(\hat\theta\) 的分布,见图中绿色部分。同样的,红色是 n=80 不断试验的分布变换。

看的出来,随着试验次数的增多 - \(\hat\theta_{10}, \hat\theta_{80}\) 都趋近于正态分布

  • \(\hat\theta_{10}\) 的分散度比 $ _{80}$ 要大,即方差要大

  • \(\hat\theta_{10}, \hat\theta_{80}\) 的均值都在 0.7

Leetcode 1029. 两地调度 (medium)

公司计划面试 2N 人。第 i 人飞往 A 市的费用为 costs[i][0],飞往 B 市的费用为 costs[i][1]。

返回将每个人都飞到某座城市的最低费用,要求每个城市都有 N 人抵达。 

示例:

输入:[[10,20],[30,200],[400,50],[30,20]] 输出:110 解释: 第一个人去 A 市,费用为 10。 第二个人去 A 市,费用为 30。 第三个人去 B 市,费用为 50。 第四个人去 B 市,费用为 20。 最低总费用为 10 + 30 + 50 + 20 = 110,每个城市都有一半的人在面试。

提示:

1 <= costs.length <= 100 costs.length 为偶数 1 <= costs[i][0], costs[i][1] <= 1000

链接:https://leetcode-cn.com/problems/two-city-scheduling

暴力枚举法

最直接的方式是暴力枚举出所有分组的可能。因为 2N 个人平均分成两组,总数为 \({2n \choose n}\),是 n 的指数级数量。在文章24 点游戏算法题的 Python 函数式实现: 学用itertools,yield,yield from 巧刷题,我们展示如何调用 Python 的 itertools包,这里,我们也用同样的方式产生 [0, 2N] 的所有集合大小为N的可能(保存在left_set_list中),再遍历找到最小值即可。当然,这种解法会TLE,只是举个例子来体会一下暴力做法。

{linenos
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
import math
from typing import List

class Solution:
def twoCitySchedCost(self, costs: List[List[int]]) -> int:
L = range(len(costs))
from itertools import combinations
left_set_list = [set(c) for c in combinations(list(L), len(L)//2)]

min_total = math.inf
for left_set in left_set_list:
cost = 0
for i in L:
is_left = 1 if i in left_set else 0
cost += costs[i][is_left]
min_total = min(min_total, cost)

return min_total

O(N) AC解法

对于组合优化问题来说,例如TSP问题(解法链接 TSP问题从DP算法到深度学习1:递归DP方法 AC AIZU TSP问题),一般都是 NP-Hard问题,意味着没有多项次复杂度的解法。但是这个问题比较特殊,它增加了一个特定条件:去城市A和城市B的人数相同,也就是我们已经知道两个分组的数量是一样的。我们仔细思考一下这个意味着什么?考虑只有四个人的小规模情况,如果让你来手动规划,你一定不会穷举出所有两两分组的可能,而是比较人与人相对的两个城市的cost差。举个例子,有如下四个人的costs

1
2
3
4
0 A:3,  B:1
1 A:99, B:100
2 A:2, B:2
3 A:3, B:3
虽然1号人去城市A(99)cost 很大,但是相比较他去B(100)来说,可以省下 100-99 = 1 的钱,这个钱比0号人去B不去A省下的钱 3-1 = 2 还要多,因此你一定会选择让1号人去A而让0号人去B。

有了这个想法,再整理一下,就会发现让某人去哪个城市和他去两个城市的cost 差 $ C_a - C_b$相关,如果这个值越大,那么他越应该去B。但是最后决定他是否去B取决于他的差在所有人中的排名,由于两组人数相等,因此差能大到排在前一半,则他就去B,在后一半就去A。

按照这个思路,很快能写出代码,代码写法有很多,下面略举一例。代码中由于用到排序,复杂度为 \(O(N \cdot \log(N))\) 。这里补充一点,理论上只需找数组中位数的值即可,最好的时间复杂度是 \(O(N)\)

代码实现上,cost_diff_list 将每个人的在原数组的index 和他的cost差组成 pair。即

1
[(0, cost_0), (1, cost_1), ... ]

这样我们可以将这个数组按照cost排序,排完序后前面N个元素属于B城市。

{linenos
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# AC
# Runtime: 36 ms, faster than 87.77% of Python3 online submissions
# Memory Usage: 14.5 MB, less than 14.84% of Python3 online
from typing import List

class Solution:
def twoCitySchedCost(self, costs: List[List[int]]) -> int:
L = range(len(costs))
cost_diff_lst = [(i, costs[i][0] - costs[i][1]) for i in L]
cost_diff_lst.sort(key=lambda x: x[1])

total_cost = 0
for c, (idx, _) in enumerate(cost_diff_lst):
is_left = 0 if c < len(L) // 2 else 1
total_cost += costs[idx][is_left]

return total_cost

转换成整数规划问题

这个问题对于略有算法经验的人来说,很类似于背包问题。它们都需要回答N个物品取或者不取,并同时最大最小化总cost。区别在它们的约束条件不一样。这道题的约束是去取(去城市A)和不取(去城市B)的数量一样。这一类问题即 integer programming,即整数规划。下面我们选取两个比较流行的优化库来展示如何调包解这道题。

首先我们先来formulate这个问题,因为需要表达两个约束条件,我们将每个人的状态分成是否去A和是否去B两个变量。

1
2
x[i-th-person][0]: boolean 表示是否去 city a
x[i-th-person][1]: boolean 表示是否去 city b

这样,问题转换成如下优化模型

\[ \begin{array}{rrclcl} \displaystyle \min_{x} & costs[i][0] \cdot x[i][0] + costs[i][1] \cdot x[i][1]\\ \textrm{s.t.} & x[i][0] + x[i][1] =1\\ &x[i][0] + x[i][1] + ... =N \\ \end{array} \]

Google OR-Tools

Google OR-Tools 是业界最好的优化库,下面为调用代码,由于直接对应于上面的数学优化问题,不做赘述。当然 Leetcode上不支持这些第三方的库,肯定也不能AC。

{linenos
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
from ortools.sat.python import cp_model

costs = [[515,563],[451,713],[537,709],[343,819],[855,779],[457,60],[650,359],[631,42]]

I = range(len(costs))

model = cp_model.CpModel()
x = []
total_cost = model.NewIntVar(0, 10000, 'total_cost')
for i in I:
t = []
for j in range(2):
t.append(model.NewBoolVar('x[%i,%i]' % (i, j)))
x.append(t)

# Constraints
# Each person must be assigned to at exact one city
[model.Add(sum(x[i][j] for j in range(2)) == 1) for i in I]
# equal number of person assigned to two cities
model.Add(sum(x[i][0] for i in I) == (len(I) // 2))

# Total cost
model.Add(total_cost == sum(x[i][0] * costs[i][0] + x[i][1] * costs[i][1] for i in I))
model.Minimize(total_cost)

solver = cp_model.CpSolver()
status = solver.Solve(model)

if status == cp_model.OPTIMAL:
print('Total min cost = %i' % solver.ObjectiveValue())
print()
for i in I:
for j in range(2):
if solver.Value(x[i][j]) == 1:
print('People ', i, ' assigned to city ', j, ' Cost = ', costs[i][j])

完整代码可以从我的github下载。

https://github.com/MyEncyclopedia/leetcode/blob/master/1029_Two_City_Scheduling/1029_ortool.py

PuLP

类似的,另一种流行 python 优化库 PuLP 的代码为

{linenos
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
import pulp

costs = [[259,770],[448,54],[926,667],[184,139],[840,118],[577,469]] # 1859


I = range(len(costs))

items=[i for i in I]
city_a = pulp.LpVariable.dicts('left', items, 0, 1, pulp.LpBinary)
city_b = pulp.LpVariable.dicts('right', items, 0, 1, pulp.LpBinary)

m = pulp.LpProblem("Two Cities", pulp.LpMinimize)

m += pulp.lpSum((costs[i][0] * city_a[i] + costs[i][1] * city_b[i]) for i in items)

# Constraints
# Each person must be assigned to at exact one city
for i in I:
m += pulp.lpSum([city_a[i] + city_b[i]]) == 1
# create a binary variable to state that a table setting is used
m += pulp.lpSum(city_a[i] for i in I) == (len(I) // 2)

m.solve()

total = 0
for i in I:
if city_a[i].value() == 1.0:
total += costs[i][0]
else:
total += costs[i][1]

print("Total cost {}".format(total))

代码地址为

https://github.com/MyEncyclopedia/leetcode/blob/master/1029_Two_City_Scheduling/1029_pulp.py

Your browser is out-of-date!

Update your browser to view this website correctly. Update my browser now

×